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ABSTRACT 

As part of the SEASAT program of NASA, a set of four 
hemispheric, atmospheric prediction models were developed 
for ECON, Inc. under contract NASW-2558. These descriptors 
applied to the four models, vyhich use a polar stereographic 
grid in the horizontal and a sigma coordinate in the vertical, 
are: 

* PECHCV - five sigma layers and a 63 x 63 horizontal grid 

* PECHFV - ten sigma layers and a 63 x 63 horizontal grid 

• PBFHCV - five sigma layers and a 187 x 187 hox'izontal grid 

• PEFHFV - ten sigma layers and a 187 x 187 horizontal grid. 

Conservation forms of the difference equations based on 
the Arakawa technique are integrated using either a fifteen 
or five minute time step on a 381 km. or 127 3an. grid (at 60° N) 
for the 63 x 63 or 187 x 187 models, respectively. Pressure 
gradient force terms are replaced by a single geopotential 
gradient on local pressure surfaces ho reduce inconsistent 
truncation error (Kurihara modification) . Stress is applied 
at the lowest level. A nonlinear pressure smoothing is used 
to help control computational noise. The horizontal boundary 
conditions are insulated slippery walls. Centered time 
differencing with time averaging of the pressure gradient 
force term in the momentum equations is used. Robert time 
filtering of the temperature and moisture solutions is used 
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for computational stability. 

The moisture and heat source/sink terms are modeled 
in a similar manner to those in the early Mints and Arakawa 
general circulation model. Terms representing evaporation 
and large-scale condensation, sensible heat exchange, 
parameterized cumulus convection cind precipitation, and 
solar and terrestrial radiation are included. Dry convective 
adjustment precludes hydrostatic instability. 

Initialization of the models is based on a pattern 
conservation technique to obtain objective analyses of the 
state parameter structure from the surface to 50 mb. 
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I. INTRODUCTION 

This document describes a set of four primitive- 
equation atmospheric forecast models for the northern 
hemisphere that were designed and developed under NASW- 
2558 for ECON, Inc. as part of the SEASAT program of NASA. 
The models are hemispheric, using a polar stereographic 
grid in the hox'izontal and a sigma coordinate in the ver- 
tical. The descriptors applied to the four models are; 

PECHGV - five sigma layers and a 63 x 63 
horizontal grid 

PECHFV •- ten sigma layers and a 63 x 63 
horizontal grid 

• PEFHCV - five sigma layers and a 187 x 187 

horizontal grid 

• PEFHFV - ten sigma layers and a 187 x 187 

horizontal grid 

The 63 X 63 and 187 x 187 grids correspond to mesh lengths 
of 381 km and 127 ]cm at 60“N, respectively. 

There are tv/o main Sections describing these forecast 
models. First, Section II, the physical-mathematical 
description, gives the governing equations and describes 
the various physical dynamical processes and their param- 
eterizatiotts that are included. Second, Section III 
describes the program structure of the various models. 

In both sections,' and particularly in Section II, the 



approach, has been to describe the five layer, 63 x 63 
model as a baseline and then to describe the necessary 
modifications to either increase the vertical resolution 
to ten layers or to increase the horizontal resolution to 
187 X 187, or both. 



II. PHYSICAL-tlATHEMATICAL DESCRIPTION OF THE 
FORECAST MODELS 

This section gives the meteorological and mathematical 
description of the set of four forecast models; in other v7ordSf 
the governing equations, numerical techniques, algorithms 
and initialization procedure V7hich constitute the theoretical 
basis of the models. 

The basic model equations and the primitive equations 
are given in Section II-A in a form applicable to a polar 
stereographic map projection v^ith an x,y coordinate system. 

Section II-B then describes the grid structure and bound- 
ary conditions for the hemispheric grid. Included is a 
description of the model levels , horizontal and vertical 
grid structure and boundary conditions. Section II-C describes 
the heating and moisture source terms which include the solar 
and terrestrial radiation and sensible heating, large-scale 
condensation and convective adjustment, both parameterized 
cumulus convection and dry convective adjustment. Section 
II-D describes the numerical procedures used to approximate 
the solution of the primitive equations. \ 

Finally, Section II-E describes the methods used to 
create an initial state for a forecast from the data provided 
by the analysis models. 
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A. 


Basic Model Equations 


The governing partial differential equations for the 
forecast models, written in flux form, are similar to sets 
used by Smagor insky et (1965) and Arakawa et ^ (1969) . 
These equations are listed below for an x-y map projection 
with map factor ra and vertical coordinate, cf. 

Momentum equation in the x direction: 


3 (tm) _ 
3t 



,3 ,uu'n-. , 3 ,uviT., 

^1T> 3? ^ 


+ n 


3 (wu) 
3a 


[11,1] 


+ TTVf 


m 


4 - 

<’'a3E + 


RT 


3ir. 


+ F + D 

X 


u 


Momentum equation in the y direction; 
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Thermodynamic energy equation: 
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Moisture conservation equation: 
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Pressure tendency equation: 
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Hydrostatic equation: 

30 _ _ RT 
3a ~ a 


In the above equations, 

u = X direction wind component 

V = y direction wind component 

TT = terrain pressure 

T = temperature 

q = vapor pressure 

w = “ ^ =* vertical velocity 

m = map factor 

f = Coriolis parameter 

F^,F^ = stress components 

D ,D - diffus-ion terms % 

u V ■■ - 

H = diabatic heating term 

Q = moisture sour-ce/sink term 

0 — gz = geopotential 
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B. Grid Structure and Boundary Conditions 

Integration of the forecast model equations is performed 
on Phillips (1957) sigma surfaces in which pressure f p, is 
normalized with the underlying terrain pressure r ir; that is: 

0=1- [II. 7] 

TT 

Figure Il-i shows the sigma surfaces used in the forecast 
models . 

The horizontal wind components;- u and v, temperature 
T and geopotentials ;r $ ;■ are carried at the levels where 
0 — 0.9, 0.7, 0.5, 0.3 and 0.1 for the five layer models, and 
at 0 = 0.95, 0.85, 0.75, 0.65, 0.55, 0.45, 0.35, 0.25, 0.15 
and 0.05 for the ten layer models.. The moisture variable, 
q, is carried at the lower three and lower six of these sigma 
surfaces for the five and ten layer models, respectively. 

The vertical velocity, yir is calculated diagnostically from 
the continuity equation for the layer interfaces (0 = 0.8 , 

0.6, 0.4, and 0.2 for the five layer models, and 0 = 0.9, 0.8, 
0.7, 0.6, 0.5, 0.4, 0.3, 0.2 and 0.1 for the ten layer models) 
and is assumed to vanish at 0 = 1. 0 and 0 = 0. All variables 
are carried at each point (i,j) of the horizontal grid. 

The hemisphere is mapped onto a polar stereographic 
projection true at 60® north for x-fhich the map factor is 
given by; 
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FIGURE . II-l : DIAGRAJ'I OF SIGMA SURFACES AND VARIABLES 
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where ^ is the north latitude. The square grid is centered 
about the Worth Pole (and extends approximately to the equa- 
tor ^ as shown in Figure II-2) . A 63 by 63 or 187 by 187 
horizontal grid system is used which gi.ves a grid distance 
of Ax = Ay = 381 km and 127 3om at 60°Nr respectively. 

The horizontal boundary conditions are similar to those 
used in the six-layer WMC forecast model described by 
Schuman et a^ (1968) . The model is bounded by rigid, imper- 
meable vertical walls placed on the grid on the next to 
outermost row of grid points, as shown in Figure II-2. 

The boundary is both insulated and slippery; that is, no 
heat or momentum exchange is permitted across the wall, 
but a parallel flow is permitted along the wall. This means 
that, along the lower boundary row where y is fixed, the 
boundary conditions can be written as: 

(iru)^ - (iru) ^; {TTv) jj^ = - (irv) ^ tn.9] 

T^ = T3? (irq)^ = (1^3)3; = ^3 

X'Tliere the subscripts refer to the first and third rows 
(row 2 being the boundary) . 




since the polar stereographic map projection becomes 
quite distorted in the corner areas of the grid and since 
these regions are not of meteorological significance for 
the forecast, the computed tendencies are truncated near 
the edge of the grid. From 15®U to the Pole, full values 
of the computed tendencies are used. Below 5®N, one-half 
of the computed tendencies are used and the amount used 
varies linearly from one-half to one between 5“N and 15“N. 
This procedure also increases the stability of the compute 
tion near the edge of the grid. 
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C. Heating and Moisture Source Terms 


The diabatic heating rates in the forecast models are 
calculated in a manner similar to that of Mintz and Arakawa, 
as outlined by Langlois and Kwok (1969) . The equation for 
the total heating can be written as 


H = «SW ^ «LW + Hp + 


ill. 10] 


where Hg^ is the heating due to absorption of shortwave 
(solar) radiation; is the heating due to longwave (ter- 


restrial) radiation; h-, is the contribution of sensible 
heat from the surface layer (cr = .9 to o = 1.0 for the five- 


layer models and o = 0.95 to a = 1.0 for the ten-layer models, 


a constant flux boundary layer) ; and is the release of 
latent heat resulting from phase changes of water substance, 
usually the latent heat of condensation. 

Section 1 describes the computation of the heating 
rates due to shortwave radiation, Hg^, longwave radiation, 
sensible heat, Hj.. Section 2 describes the large- 
scale condensation process which contributes to H^. Section 
3 describes the moist and dry convective adjustment proc- 
esses which also contribute to 
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1 . Solar and Terrestrial ]^diatlon and Sensible Heating 


Diagnostic calculations of the heating; i.e*^ the tem*^ 
perature tendency, AT/At, represented by Hg^r and Hj., 
are computed every hour and saved for use over the next 


one-hour forecast period. The temperature tendencies are 
of the form 


AT _ g AF 

2^ “ Cl Sp 

p 


[ 11 . 11 ] 


where g is the acceleration of gravity; is the specific 
heat of air at constant pressure and ^ is the flux diver- 
gence. Flux computations depend on the vertical distribu- 
tion of temperature, moisture and pressure above each sur- 
face grid point; i*e., only vertical fluxes are computed. 

The computations of Hg^, and Hj, are described in 
Subsections 2, 3 and 4, respectively. 
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a. Solar Kadiation 


The incoming solar radiation at the top of the model 
atmosphere is given by 


S - Sq COS y [11.12] 

2 . 

where Sq is the solar constant (taken to be 1.92 eal/cra min.) 
and y is the zenith angle which is a function of time of day. 
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latitude and season. 

The zenith angle Is given by 

cos y “ -sin4> sin 6 + cos(|) cos 6 cosh [II, 13] 

where 41 - latitude (0® at Equator, 90® at North Pole) 
h = hour angle 

6 = declination angle of sun 

and 

sin5 = 0.39785 sin [4. 88578+0. 0172*DAY 

+ 0.03342 sin (0. 0172*DAY)^0. 001388 cos (0.0l72*DAY) 

+ 0.000348 sin (0,0344*DAY) 

- 0.000028 cos (0.0344*bAY)] 

where DAY is the nuaanber of days measured from the first of 
the year. 

In order to compute the depletion of the solar beam, S 
is initially partitioned into two segments. The radiation 
with wavelength shorter than .9 microns {65.1% of the energy) 
is subject to scattering by the atmosphere, but not to absorp-» 
tion. Wavelengths greater than .9 microns are subject only to 
absorption. The calculations that are outlined below follov? 
the scheme constructed by Joseph (1966) . 

In order to compute the radiative fluxes, a simple form 
of middle cloud is assumed to exist between a = 0.4 and a~ 0.8 
An empirical formula for middle clouds: 

P 

CL = 1.3 (|-)- - [0.46 + 0.6 (-%r-^)] [1 - (1-*)^] [II. 14] 

®s ^ ^2 ®s ^ 
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gives the amount of cloudiness/ CL in terms of the relative 
humidity and pressure at u = C. C is taken to be 0.7 for 
the five layer mciels and 0.65 for the ten layer models. 

= 400 and Cg = 300 for the five layer models, and 300 

and G 2 ~ 350 for the ten layer models. Limits of 1.0 and 0.0 
are set on the range of values given by equation [11.14]. 

The absorption of incoming solar radiation for the upper 
layer (a = 0 to d = 0.6) is given by: 


A 2 = 0.349St{l~CL)Fj^ + CL(l-a^)F2l 


[Ii.15] 


and for the lower layer (a = 0.6 to o = 1.0) by: 


*= 0.349S[(1-CI.)F3 + CL(l-a^)F^] 


where 


F, = 


0.271 [ESSm] 
*^cosu ■' 


(U), 0.303 


F» =0.271 [2.91 


0.303 


F^ = 


F, = 


jEga, 0.303 . 

1 [2.91 pw(T)]®*^°^ 


The symbol pw is the precipi table water given by: 


1 t^2 

P« = S /„ 0-622 


1000 


dp 


[11.16] 


[11.17] 

[11.18] 

[11.19] 

[ 11 . 20 ] 


[II. 21] 


and U refers to the amount of preeipitable water above <s = 0.6 
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u 


G 


f 1 


and T refers to the total amoxmt of precipitable water , 

These absorption functions are from Manabe and MoXler |1961) . 

The clouds' albedo, a_, is set at 0.5. Finally, the amount 

'C 

of solar radiation that reaches the surface , , is needed 

for the calculation of sensible heating over the land areas. 
Following Hints and Arakawat 

0.651S(l-a ) 

S„ = a-a) ( tl-CL) r0.349S - A, - A. + . 
y y s g 

0.349Sl(l-aJ - (1-a )(F. +F.)] 1-a 

+ CL { a + 1=^ 0,651S}) 

e g sq g 


[II .22] 


O 


where a , the ground surface albedo, is given as a function 

g 

of cost over the sea and as a function of snow cover over 
the land; a__ is the albedo of the cloudy sky and a is the 

SO ^ 

albedo of the clear sky and is given by: 

P_ 


ttg = 0.085 - 0.10727 ln[siny (j-^) ] 


[11.23] 


b . Terrestrial Radiation 


i ('■. 


The infrared fluxes at 0 - 1,0, 0.6 and- 0 . 2 are computed 
in the manner suggested by Mints and Arakawa and are given by: 


O 


1.0 


0.6 


*■ 1.0 ♦ 


^ 0.6 


0»8(1-CL) - Fi q) 

1 + 1.7S[pw(L)l°*^^® 


[11.25] 
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and 



^ 0.2 " %.2 


o.8(i-ei.) 


1 + 1.75Ipm(I.> + 1.2 pw(U)l 


KtlT 


where St is the Stefai^~Boltsmanh constant and U L to 

tlte amount of precipitable water above and below a » 0 . 6 / 
respectively. 

lit it it 

Pi QP Fq 6 ^0 2 fluxes computed by assuming 

the surface temperature is given by a downward extrapolated 
temperature j C3 Tq g “ Tq 7 ) for the five layer models 

and Tg = ^ (3 Tq gg “ ^0 85^ layer models. The 

factor 0.8 is an empirical correction to account for the 
effect of CO 2 . 

In equation [11.24]/ T„ is given over the oceans by the 

■ ■ • y . 

sea surface analysis/ and over the land areas it is given by 
the solution of a heat balance equation to be explained in 
Section II-C-l-c on sensible heat transfer and evaporation, 

it 

To compute the approximate fluxes / F , the method devel*' 
oped by Danard (1969) is used with an emissivity function 
for water vapor from Kiahn (1963) . Thus, 

= St ( (1-CL) ( -T^+Tg^Q$ (Xl)+T |^4 lt(X2) 

CUtJ grT|-yg^g)*(3(4)r ^ 11^^21] 

pj^g » St{l-Cb) t^+tT|-Tj^glt|3^^^ [11.28] 
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c. 


Sensible Heating and Evaporation in the 
Arakawa-Type Planetary Boundary Layer 


Sensible heating is computed as a function of the dif- 
ference of the surface temperature, and the temperature 
of the air near the surface, T . Similarly, evaporation is 
computed as a function of the saturation specific humidity 
at the surface, Q , and the specific humidity near the sur- 

s 

face, Q , except over land, where a Bowen ratio is computed 
to determine the flux of moisture. Difficulty arises in 
the forecast models since near-surface temperature, wind and 
moisture are not predicted. To surmount this difficulty, the 
technique of Mints and Arakawa is used. 

The air near the surface; i.e., a = 1.0 to 0 ~ 0.9, 


in the five layer models and a = 1,0 to 0 = 0.95 in the ten 
layer models, is assumed to be a very thin boundary layer 
which does not absorb or store a significant amount of heat 
or moisture. The fluxes of sensible heat and moisture are. 


therefore, the same through the top and the bottom of this 
layer. Generally, for heat transfer, an equation of the 


fojnn: 


r = -K 


49 

H AZ 


III. 38] 


is used where 0 is the potential temperature. Equation 
[11.38] is modified by Arakawa to 

* 

" = Pl.O c 


K 


'p a*t0BL"®S* 


[y - j 

^’c Az ^ 


III. 39] 


Az 
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where AZ « ^ . In equation i ll * 39 j > the expression 

for Kg is allowed to vary with the stability of the boundary 
layef . K , a and y_ are empirical constants. The eonstant 
Vw allows for a somewhat more realistic variation of r between 
day and night . The subscript s refers to the ^iue of a 
parameter extrapolated to the surface from the two lowest 
sigma levels. 

The sensible heat flux through the bottom of the 


layer is given 




where#, for this computationf 


- 0.8 -hv|j,] V2 + BO 


where Gu, the gustiness factor# is 2,2 m/sec, A drag 


coefficient ^ 0v002 is used. - ^ 

oyer the oceans, T„ can easily -be obtained fromv 
eguations a^^^^ # Since the sea 


erathre T^, is specifiedv Knowing T^j,, a heaMng rhid 

.'tb-."seh|^bieCheai ' eigm'a-' . layer;;, 


■'face ■enefgyi'iaaianca-- 




'Slrl?.' 








in which only the heat storage has been neglected. Data 
from Budyko {given in Sellers 1965) has been used to deter- 
mine r, a Bowen ratio, as: 

r = 9.6 (sin({))^ - 7.93 (sinc}>) + 2.0 [11.43] 

in equation [11.42], the amount of solar radiation reaching 
the ground, S^, is given by equation [11.22] , and the infra- 
red flux from the surface, q, is given by equation [11.24]. 
Equation [11.42] is modified over ice-covered ocean to: 

(l+r)Hj. “ Sg + Q = B(T^ - Tg) [11.44] 

2 

where B = 697.83 erg/cm sec °K is the ice Gonduction coef- 
ficient for three-meter thick ice and T, = 271. 2®K. 

Evaporation over the ice-free oceans is determined in a 
similar manner to the heat flux. The equation resulting is: 


where 


^1.0 S^s , V 

E = (qg - q^) 


1 + 




Az 


= 


K 




1 + 


Az 


q 


s 


0.622 e 


P_ - e. 


[11.45] 


[11.46] 


[11.47] 
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[11.48] 


Sg = exp(Ae - 

Ae = 21.656 and Be - 5418. 

The details of the derivation of equation [11.45] are given 
by Langlois and Kwok (1969), 

Using the Bowen ratio from equation [11.43] allows the 
determination of the evaporation for land areas and ice- 
covered ocean as; 


E = M/r [II,49j 

The expression for r gjven by equation [11,43] has the effect 
of making the land surface at low latitudes a greater source 
of moisture than those at higher latitudes. 

2. Large-Scale Condensation 

The large-scale condensation mechanism operates under 
the premise that synoptic scale motions do not contain super- 
saturated states. If a state of supersaturation is discovere 
(which may be due to excessive evaporation from the underlyir. 
surface or due to the effects of horizontal or vertical advec 
tion of moisture) , then sufficient moisture is removed (as 
precipitation with the latent heat of condensation) to achiev 
a saturated state at a warmer temperature. The technique 
suggested by Mintz and Arakawa as outlined by Langlois and 
Kwok (1969) is followed in the forecast models. 
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A solution to tho equation 


F(T) == Cp(T-T^) - Llg^-qg(T)] =0 [11.50) 

is sought; v/here is the suf>ersaturation mixing ratio 
associated with the temperature T cuid g is the saturation 
mixing ratio associated with the new temperature T. L is 
the latent heat of condensation which is taken to be 600 


cal/gm . 

Equation [11.50] is a non-linear equation. A second- 
order Hev/ton-Raphson approximation to the solution is used, 


which is given by: 


T 


T 


o 


F(T^) 

F''"(T^) 


[1 + 


F' * (T ) 

o 

2 


F(T ) 
SL— } 

F' (T^)^ 


[11.51] 


Expressions for the derivatives of P(T^) must now be 
generated. F'{T) may be written as: 


F-TO = 1 + ^ q^(T) 

P 


[11.52] 


where 


ee 


q (T) = — 2- 

^s'‘ p-e. 


and e is given by equation [11.48]. Differentiating 

q (T) and ignoring dp/dT gives; 
s 


F* (T^) 


= 1 + 


L sp Be 
C~~ ■> 2 

p V 


(T ) 

s o 


[11.53] 
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which may be further simplified to give 


F' (T ) = 1 + ^ Y 


[11.54] 


where : 


Y 


Beq^(T^) 



[1 + 


“Js '’’o> ' 


[11.55] 


The second derivative F*'(T^) is given by: 


F' ' (T ) 
o 


L dY 

C dT 
P 


[11.56] 


which, after coiisiderable algebraic manipulation, can be 
written as : 

F' • (T^) = (| + I q^) - 1] [11.57] 

p o o 

Using equation [11.51], the change in temperature may 
he obtained and from this, the change in mixing ratio can 
be computed by: 

Aq “ AT [11.58] 

The air is warmed by the amount AT and moisture is removed 
by the amount Aq from the supersaturated layer. 

Since moisture is forecast at the three lower levels 
in the five layer models and at the lower six levels in the 
ten layer models, a precipitation algorithm is required. 
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This algorithm proceeds from the lowest level, to the highest 
level, and is outlined below, 

• Examine the layer for super saturation. If we are 
at the lowest sigma level and the layer is supersaturated, 
then precipitation occurs? otherwise, proceed to the next step. 

• Examine the next lower layer for saturation. If 
the layer is saturated, then we allow the precipitation to 
-pass through the layer to the surface where it is added to 
the total at that particular grid point; otherwise, proceed 
to the next step. 

• If the next lower layer was initially unsaturated and 
the addition of the precipitation from the layer above results 
in supersaturation, we repeat the calculations summarized 

in equations [11.51] and [11.58]. The moisture removed is 
treated as precipitation; otherwise, proceed to the next step. 

• If the additional precipitation does not result in 
supersaturation at the next lower level, we add the moisture 
to the layer and decrease the temperature accordingly. 

Finally, precipitation is given by the equation: 

Pr = In(cr) Aq (cm) [11.59] 

where r = 1.0/0, 8 or 0.8/0. 6 for the levels in the five layer 
models, and r - 1.0/0, 9, 0.9/0. 8, 0. 8/0,7, 0.7/0. 6, or 0.6/0. 5 
for the levels in the ten layer models. 
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3 . Convective Adjustment 


In the atmosphere, the direct consequence of vertical 
static instability is the initiation of small-scale vertical 


convection which is very active in transferring heat and 


moisture from the surface through the lowest layers of 
the atmosphere. This process is very instrumental in main- 
taining the overall potential energy of the atmosphere against 
its destruction by radiative absorption and emission. 

In the forecast models, the primary purpose of the con- 
vective adjustment mechanism is the removal of vertical in- 
stabilities, whether they be generated by a large-scale dy- 
namical process or by heat and moisture transfer from the 
surface to the lowest computational level of the model. If 
the vertical instabilities are allov/ed to grow, they will 
promote vertical convection whose scale is dictated only by 
the grid length in the model. This will lead to computational 
instability and destroy the solution. In order to control 
this type of instability and also to simulate vertical con- 
vection in the atmosphere, a convective adjustment process 
consisting of two parts, a parameterized cumulus convection 
and a dry convective adjustment, is performed at each time 
step. 
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Parameterized Cumulus Convection 


a. 


Following a method devised by Arakawa et al (1969) , 
three types of parameterized cumulus convections are con- 
sidered. The first type is that occurring between levels 
2 and 3; the second type is that occurring between levels 

1, 2 and 3; and the third type is that occurring between 
levels 1 and 2. Precipitation may result from types 1 and 

2, but not from type 3. Figure II-3 schematically illustrates 
these three types of parameterized convection. 

Both the five layer forecast models and the ten 
layer forecast models use the same parameterized cumulus 
convection levels as shown in Figure II-3, In the case of 
the five layer models, the convection levels coincide with 
sigma surfaces, while in the case of the ten layer models, 
two adjacent sigma levels are combined to form a composite 
convection level with an averaged temperature and mixing 
ratio. Should cumulus convection then take place, the 
changes in moisture and temperature are apportioned between 
the two model layers on the basis of their mixing ratios. 

Three parameters used to indicate which type of 
convection can occur are the energy integrals 

S = C T + gz [11.60] 

P 

h =* Cp T + gz + Lq = S + Lq [11.61] 

E = Cp T + gz + Lqg = S + Lq^ [11.62] 
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E is constant along a moist pseudo- ad iabat and S is constant 
along a dry adiabat. Thus, comparison of h and E for var- 
ious levels allows for a sufficient test of conditional 
instability. 

Type 1 convection will occur if 


*^2 ’ 


If E 2 = E^/ a moist adiabatic lapse rate would exist ^ 
between levels 2 and 3. Since h 2 < E 2 , inequality [11.63] 
implies that a conditionally unstable lapse rate exists. 
Furthermore, inequality [11.63] can be rewritten as: 


[11,63] 


*52 ^ L ^^3 " ^2^ 


Type 2 convection will occur if 


E^ > h 2 and h^^ > E^ 


[11.64] 


[11.65] 


Conditional instability does not exist here between levels 
2 and 3, but does exist between levels 1 and 3. 

Finally, type 3 convection will occur if; 


E^ > h-i > E2 


[ 11 . 66 ] 


Levels 1 and 2 are conditionally unstable; levels 2 and 3 
are not. 

The procedure now is to write do^vn equations describing 
the budget for the moisture, q, and the dry static energy, S, 
for each type of convection. Before this can be done, the 
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energy budget for the lower part of the cloud is considered, 
neglecting any accumulated storage of energy within the 
cloud. 

We use type 2 convection as an example. When the 
entrainment coefficient, ri # is greater than 1, entrainment 
exists and the parameter E in the clouds, E , is given by 
a mixing of hj^ and in the lower part of the clouds, as 
follows: 

' "2 n ”'1 • ‘'2> 

because h is approximately conserved with respect to an 
air parcel. At level 3, 


= Vc3 ^ 9^3 ^ 


(11.68] 


whereas, for the environmental air: 


Therefore , 


Ej . CpT, = 


3 = gZj + Lq^CTj) 


1 

T - T = ^ ^ — 

"c3 ^3 1+Y„ C 


[ 11 . 69 ] 


[ 11 . 70 ] 


where 


= L ^^s 

*^^3 C 9T 


[ 11 . 71 ] 


For type 1 and 3 convection, is given by and h^^, respec- 
tively. From equation [11,70], one can obtain approximately; 


^3 


*c3 ^s3 1+Y- 


[11.72] 


For types 1 and 2, the convective terms in the moisture and 
dry static energy budget can be written in flux form as: 


- r fn + “^3 _ ‘32+^3. 

g 3t L 

Ap, ^2 1 

^ = >'= ‘- V ^ ^ I 


Ap, 9q, gp-gi 

at^ = C(-^^) (type 2 only) 


= 0-=3 , " 3 -" 2 , 


g 

at 

AP3 

883 

g 

at 


883 

g 

8t 

Ap^ 

88^ 

g 

3t 

3 convec 

Ap2 

aqj 

g 

8t 


^ 3"®2 1 ^ 2~^1 

nc [-i^ + i (^^)) 


^2-^1 


= C ( ■ 2 ^ (type 2 only) 


g 9t 


q]-q2 

c (- V ") 
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g at 





[11.80] 


where ti is the entrainment parameter and C is the convective 
mass flux. 

Following Arakav/a (1969) we now state that the 
solution of the above differential equations has a charac- 
teristic time scale, t, such that: 


It 


(h-E) 

T 


For type 1 convection, we can then obtain 


[11.81] 


gtj^ ‘2+Y ' 




,1+Y. 


[11.82] 


3 I(h^-S3)-f( "^'3 ) (S 3 -S 2 )] 


where is given by the empirical relations; 

- 3600 [l+0.667(h2-E3)] ^ [11.83] 


and 


= 1200 [22(h2-E3) - 7] ^2~^3 - 


[11.84] 


For type 2 convection, a value of 1®K was used for T - 

C ^ t 

therefore, we obtain: 


n = 


^ l “^2 


E3+Cp(l+Y3)-h2 


and 


gi2 


(hi-E ) 


S3—S2 ^1~^2 

[n(l+Y 3 ) (C +-^^) + (-~-^)] 


[11.85] 


[ 11 . 86 ] 
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where T 2 is given be the empirical relations 

T 2 = 3600 Il+2.25(hj^-E3)] < 0.5 [11.873 

and 


T 2 = 3600 [15.75(h^-E3)-5.75] > 0.5 [11.88] 


For type 3 convection, one can obtain: 

1 


G = - 




qi-q 


1 '=‘2, 


[ (~2^> (S;j^-S2) “L (— 2 • ^ ^ 
where is given by the empirical relations: 


and 


[11.39] 


T 3 = 3600 [1+0.667 (hj^-E 2 )] [11.90] 


T 3 = 1200 [22(h^-E2) - 7] [11.91] 


Finally, precipitation from types 1 and 2 convection is 
given by 


- t — -T eS rr -^ ’ 


b. Dry Convective Adjustment 


The dry convective adjustment process consists of 
checking the model atmosphere for static stability and mo- 
difying the temperatures when necessary to achieve static 
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stability. This is accomplished by first computing the 
vertical profile of potential temperature, 0, by: 


0 = T <100_Qt 

' Tr0,J 


[11.93] 


where tt is the terrain pressure, is the k sigma level 
and < = R/Cp 0.2858. Starting at the bottom of the profile, 
we then test to see if 0j^ > + 0g» where 0^ = 3®K. If 

this is the case, the layer is stable and no temperature 
modification is necessary. If it is not the case, then the 
potential temperatures from layer 1 through k are modified 
to achieve a stable profile, using the average potential 
temperature and a fraction of the temperature difference 
0- per layer. (The fraction depends on the number of layers 

Of 

whose temperatures are being modified.) Th-s process is 
continued until the entire profile has been checked for 
unconditional stability and modified, if necessary. The 
revised temperature profile is then recovered using equation 
[11.93] . 
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D. Numerical Techniques 


) 


The numerical techniques used to approximate the solution 
of the forecast model equations are given in this section in 
two parts. First, the spatial finite differencing tech- 
niques are described. Then a description of the time dif- 
ferencing technique is given which includes an explanation 
of the pressure-gradient force averaging technique , Robert 
time filtering and the non-linear pressure smoother. 


1. Spatial Finite Differencing 


Energy conserving (in space) difference equations 
based on the Arakawa technique (1966) are used to approxi- 
mate the spatial derivatives. Thus, the finite difference 
equation for the momentum equation in the x direction is 
given by: 

-V [-I.(u)* + f(t*v*) - (ff);]"" - [11.94] 


3 (iru) 
3t 


and in the y direction by; 


3 (irv) 

at 


n 

'V. 

* 


[-L(v)* 


f (Tr*u*) 


(M) ] 


n 


Fy"^ + ^ [11.95] 


where the subscript * denotes the point (i,j,k)? the super- 
scripts denote the time level; L is the Arakawa advective 
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operator; f is the Coriolis parameter; the A notation 
refers to the Kurihara pressure gradient; and the bar 
denotes time averaging as explained in the next section. 


F and F are the frictional terms ; D and D are the 

X V 

lateral diffusion terms. 

The Arakawa advective operator is given by; 

2 




IT . . 

Jii2 


"** 2Aa ^\+l''‘'*'k^ “ '*'k~l^'^k'^'^k-l^ ^i, j 


where a = 


& = 


UTT 

m 

vir 

m 


m = map factor 
d = grid mesh size 


The pressure gradient terms ^ and ^ are computed 
using a type of modification suggested by Kurihara (1968) . 
Instead of evaluating both the gradient of the geopotential 


distribution on a sigma surface and the terrain pressure 
gradient (invariant with height) as would be required by | 

equations [II. 1] and [11.2]^ the value of the geopotential 
is interpolated at neighboring grid point locations to the 


same pressure surface as at the computation point, allowing 
computation of the geopotential gradients on pressure surfaces 


The interpolation formula used is a form of : 




dim 


-RT 


[11.97] 


This procedure cf differencing $ on local pressure surfaces 
reduces the inconsistent truncation error over high, irreg- 
ular terrain. 

Frictional dissipation is accounted for in the gross 
boundary layer (u = 1.0 - 0.8 for the five layer models and 
a = 1.0 - 0.9 for the ten layer models) only. The stress is 
assumed to vanish at the top of this layer. The stress 
terms are the same as those given by Shuman and Hovermale 
(1968) , except that the surface wind speed is obtained by 
extrapolation from the boundary layer sigma surface. Thus, 
the friction in the x direction is given by: 


,n-l 


!fevn-i 

AuR s 


(S) 


n-1 


BL 


[11.98] 


and in the y direction by 


n-1 _ „n-l /irVv 
AaR ^ 


[11.99] 


where = 0.3[u 


BL 


y2 ]V2 
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and Cjj - 0.0015 for sea level and 0.0025 for non-sea level 


terrain 


N. 


The lateral diffusion terms are given by 


2 

m, . 


D = D 

^2 

D„ = D (ttV^v)*"^ 


[II. 101] 


[ 11 , 102 ] 


where the Laplacian is 


= Y.., . + Y. . . + 1?. . + y. . . -4H'. . [11.103] 

'* x+1,3 x-1/3 1/3+1 1/3-1 1/3 


and the diffusion coefficient, D, varies with latitude according 

6 2 

to the relation lx 10 (1-sin 4>) • 

The change in the terrain pressure is computed from 


an approximate form of the continuity equation as: 


[11.104] 


K 2 

W ]^1 2d ^“i+1, j ,k“°‘i“l, j ,k'''*^i, j+l,k”^i, j-l,k^ 


The thermodynamic energy equation is approximated by: 


-^4?^ I - L(T)* 


[11.105] 


+ ( 

cp 


^''k'^'^k-l^i, j 


™i i 

2d” '^^i, j /k^^i+1, j ,k“^i-l, j ,k^ 


^ ^i/ j /k^’^i, j+l,k“^i,j-l,k^ ^ ^ 


^ i,j 
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Vertical motioas for non-material surfaces are obtained 
from a form of the continuity equation: 


Aa / 

/ j 


tll. 106 ] 


m . 




The hydrostatic equation is given at the lower level by: 




[11.107] 


where T* is extrapolated from the lowest two levels by 
1.25T^ - O. 25 T 2 and is the terrain geopotential. For 
subsequent levels 


$ = $ _ BiB. (T + T ) 

^k+1 ^k 2 ^ 0 ,. ^^k ^ ^k+l' 


[II.1D8] 


Finally, the moisture equation is approximated by 


9 (frq) 
at 


n 


'V - L(q)J + Q*TTi"j 


[11.109] 


2. Temporal Finite Differencing 

A centered, or leap-frog, difference scheme is used 
to approximate the time derivative in the forecast model 
equations; that is 


3A 

at 


n 


2At ^ 


[II. 110] 
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In order to lengthen the permissible time step used in 
the computation, a numerical technigue known as pressure 
gradient force time averaging is used in conjunction with 
the centered time differencing. This procedure was suggested 
by Shuman (1971) and has been successfully used at the 
National Meteorological Center. Pressure gradient force 
time averaging requires arranging the integration of the 
primitive equations into two parts. First, the pressure 
tendency, thermodynamic energy and moisture equations, in- 
cluding the diabatic terms, are integrated forward to the 
(ntl) time level. Then the momentum equations are integrated 
forward in time to the (n+1) time level, but the pressure 
gradient force terms Ctefer to Equations [11.94] and [11.95]) 
are time averaged according to the relation; 

= (l-a)FGF’^ + + PGf”'*'^) [II. Ill] 

For stability, 0<a<0.5, and it can be shown that when ci=0.5, 
a time step twice as long as when a=0 can bo used. However, 
a weak eomputationa’’ ’ ''c _ajjility is present and a value of 
slightly less than one-half should be used. The value 
a=0.49 is used in the forecast models. 

A centered difference scheme is neutrally stable in 
that there is no damping, only phase shift, of the various 
frequency components present in the solution of the difference 
equations. In order to achieve damping of the shorter wave 
length components to avoid non-linear instability and, in 


particular, to damp the computational mode, the simple leap- 
frog scheme is modified slightly to time filter the solutions 
of the thermodynamic and moisture equations. This procedure 
is knovm as time filtering and was originally used by Robert 
(1366) . 

Assume that the partial differential equation that we 
wish to approximate is 

= F(A) [11.112] 


Then, using the leap-frog difference scheme, we have 

= A^“^ + 2 At F(A) [11.113] 

This step is modified by Robert to become 

(A')^"^^ = a”"^ + 2At F{A,A') [11.114] 

= d-a){A’)^ + I + (A')”'^^] 0<<^<! [11.115] 


where we time-filter the solution. It can be shown that this 
procedure adds damping to the solution and heavily damps 
the computational mode that is present in the solution of 
the finite difference equations. 


If oAt<l is the linear stability criteria for the 
simple leap-frog difference scheme, then the addition of 


Robert time filtering will change the stability criteria to; 
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III. 116] 


(oAt < [ 


-g 

2+g 


1/2 


A value of g=0.1 is used to time filter the solution of the 
thermodynaniic and moisture equations. Mo time filtering of 
the solutions of the momentum equations is done sin.:e lateral 
diffusion terms are included in these equations. 

The terrain pressure tendency equation [II. 104] is 
essentially an ordinary differential equation; thus, rather 
than using Robert time filtering of the solution to control 
small-scale noise and solution separation, a non-linear 
smoothing operator is applied every time step. This smoothing 
operator has been described by Oliger and Wellck (1970) and 
is given by; 






- |p§i^ - - P^i,j-i)l Iiiai7] 
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where ifiPsI is the maximum absolute difference between 

Ps at the point (i,j) and Ps at the four immediately neigh- 
boring points and e is the smoothing coefficient which has 
been determined to be of the order of 0.04 for control of 
computational noise. 

The smoothed surface pressure is then given by: 

(smoothed) = pf^^ + APs [11.118] 

The terrain pressure is reduced to sea level pressure and 
vice versa through the relations 

pE = u** + Att'^ [11.119] 

and 

Att’^ = Air. [1+(J^ - 1) 

^ "o 

where ir^ is the initial terrain pressure and Air^ is the 
initial difference between sea level and terrain pressure. 
These relations are also used to create surface pressure 
fields for output. 

Due to the interaction of these various computational 
devices described above (pressure gradient force time av- 
eraging, Robert time filtering, lateral diffusion and non- 
linear pressure smoothing) a doubling of the permissible 
time step is not possible. However, a fifteen minute time 


,ir^-700. , 
^ 300 ^ ^ 


[ 11 . 120 ] 
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step for the 63 by 63 models as opposed to a ten minute 
time step and a five minute time step for the 187 by 187 
models is possible which results in a considerable saving 
of computer time over simple centered time differencing. 

Finally, the integrations are restarted after each 
output cycle (normally every twelve model hours) using an 
Euler-backward or Matsuno time step. This integration 
step is given by: 

^ F[A*^1 [11.121] 

and 

+ At F[ (AMn"*"!] [11,122] 


E. Initialization 

In order to be able to forecast the meteorological 
conditions at some time in the future using a primitive 
equation model, the conditions at the present time must 
be known; that is, an initial state must be specified for 
the variables forecast by the PEM. This initialization 
process is accomplished by using analyzed data produced by 
a set of analysis programs and synthesizing this data into 
a complete specification of the initial state for the PEM, 
The following sections describe the method of inter- 
polating the analyzed data fields to the forecast model 
sigma surfaces. The specification of various other required 
data fields is also described. 


1. Constant Fields 


Constant data fields are of two types — those that 
are truly constants (sin of latitude, map factor, tendency 
truncation coefficient and terrain height) and those that, 
while not constants, are held constant during the course of 
a forecast (sea surface temperature, land-sea-ice indicator 
and albedo) . 

For a polar s ^ ideographic grid, true at 60®N and centered 
about the North Pole, the sin of latitude is given by: 

sxn (j) - gyp + gT III. 123] 

where 

STP = (31.205)^ for the 63 by 63 grid 
STP = (93.615)^ for the 187 by 187 grid 

and 

ST = (I-POL)^ + (J~POL)^ 

POL =32 for the 63 by 63 grid 
POL = 94 for the 187 by 187 grid 


The map factor is given by 

- 1 t sin 60** 

^ 1 + sin (|» 
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The tendency truncation coefficient used to reduce the 


tendencies of the variables in the outer regions of the 
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grid is given by; 


re = 1.0 for <|) > 15° 

= 0.5 for (|} < 5° 

and varies linearly between one-half and one for 
5“ and 15°. 


[11.125] 


between 


The terrain height is derived from Scripps Institute of 
Oceanography data and has been area averaged and gradient 
limited. That is, the terrain gradient at any grid point is 
constrained to be less than 2000 m per 381 km, the basic 
63 by 63 grid size. 

The land-sea- ice indicator is also derived from the 
geography and varies seasonally as ice masses expand and 
contract . 


The albedo, ALB, of the earth's surface is obtained 
from the seasonal albedo charts of the northern hemisphere 
developed by Posey and Clapp. The open sea albedo is mod- 
ified in the heating computation as a function of the zenith 
angle of the sun, vi/ according to the relation; 

ALB' = ALB + 0.54 (0.7-cosu) [11.126] 

The sea surface temperature is one of the analyzed data 
fields produced by the pattern conservation technique analysis 


programs . 


2. Terrain Pressure 






The terrain pressure, n, is determined from the ana- 
lyzed sea level pressure, P , the terrain height, Z and the 

S Jl 

heights of the two analyzed pressure surfaces bounding the 
terrain height. That is, if the terrain height lies be- 
tween Zt and Z„, which correspond to the heights of the pres 
sure surfaces, and P^, then the terrain pressure is given 


by 


Pt X Z -Z 

ir = P., where X = ^ [11.127] 


The difference betv;een the sea level pressure, P , and 
the terrain pressure, -n , is also computed and stored for use 
in reducing the forecast terrain pressure to give a forecast 
sea level pressure. 

3 . Temperature 


The sigma surface temperatures are determined by loga- 
rithmic interpolation between the analyzed pressure surfaces. 
The terrain pressure, tt, and the temperatures of the analyzed 
pressure surfaces, are used to determine the temperatures of 
the sigma surface from the follov/ing relation; 




11-44 


That is, the pressvires, and P^, are found which bound 
the sigma surface, air , and along with the corresponding 
temperatures, and T„, are used in equation [11.128]. 

Li U 

4. Geopotentials 


Having computed a set of sigma surface temperatures, and 
having obtained the terrain geopotential, the hydrostatic 
equation can be used to compute an initial set of geopoten- 
tials. That is, equations [II. 107] and [11.108] are used to 
compute the sigma surface geopotentials. 


5 . Winds 

The initial winds are obtained from an analysis pro- 
cedure which gives analyzed winds on pressure surfaces. 

These winds are given on seven pressure surfaces: 900, 700, 

500, 300, 250, 150 and 100 mb and are then vertically inter- 
polated linearly in the logarithm of pressure to the sigma 
surfaces using equation [11.128] . 

6 . Moisture 

An important part of the initialization process is the 
specif icatipn of the initial moisture distribution. The 
forecast models carry moisture at the lower three levels in 
the five layer models and at the lower six levels in the ten 
layer models and would, thus, require a large amount of ana- 
lyzed moisture data for initialization, Sinee moisture data 


are virtually non-existent over the oceans, a parameterized 
approach to the specification of the initial moisture dis- 
tribution has been taken. 

Studies by Kesel and Levrit (1974) showed that existing 
analysis/procedures for moisture initialization do not nec- 
essarily result in an optimum moisture distribution. It 
was further shown that improved forecasts could result from 
moisture initialization based upon the relative vorticity 
distribution at each model level. 

ODSI's studies have shown that the initialization of 
moisture has a significant effect on deepening lows, but 
has little effect on filling lows or highs. It was shown 
that by basing the initial relative humidity on the inten- 
sity of the geostrophic relative vorticity, specific regions 
of properly posed relative humidity v=5lues will more likely 
coincide with actively developing systems. 

Therefore, a procedure based on the intensity of the 
geostrophic relative vorticity is used to define the initial 
moisture distributions on the 1,000, 900, 700, 500 and 300 
mb pressure surfaces. These values are then interpolated 
linearly in the logarithm of pressure using equation [11.128] 
to the sigma surfaces. This procedure for each of the four 
pressure surfaces can be outlined as follows: 
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a. Compute the geostrophic relative vorticity 
for each point in the field using the relation 

RV = v^Z = i~| + ~ [11.129] 

■bY 

where 2 is the height in meters of the pressure surface. 

b. Compute the mean and standard deviation of 
the relative vorticity for the field using the relations 


1 N 

M = i 1 RV. . 
N 1,3 


[11.130] 


SX) = [i X 


[11.131] 


c. Normalize the relative vorticity field according 
to the relation 


RV. .-M 
RVN. . - 


[11.132] 


d. Compute the relative humidity for each point 
based on the normalized relative vorticity values using the 
empirical relation 


^if j 


=40+30 




[11.133] 


with the upper limit of relative humidity set at 90% and 
the lower limit set at 10%. 
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III. 


COMI>UTER PROGRAM DESCRIPTION OF THE 
FORECAST MODELS 


The designer of any large, time-consuming computer 
code is faced with numerous problems in trying to design an 
efficient and readable program. First, there are the prob- 
lem constraints on how many grid points are desired, how 
many variables are to be carried at each point, exactly 
what must be computed, what are the inputs and what are the 
desired outputs. Second, there are the hardware constraints 
— the size of main memory, the type of secondary storage 
available (rotating or extended memory) , the number and 
characteristics of the I/O channels and the speed of the CPU 
Third, there are the software system constraints — exactly 
what will the system allow one to do and what sort of high 
level language is available and what are its characteristics 
Finally, there is the requirement that the resulting code 
should be efficient, readable and easily modified, both from 
the standpoint of computational changes and portability 
between computers. 

What follows is a description of the realization of 
these various design objectives, including the various com- 
promises that were required. There are four forecast models 
described that range in data base size from approximately 
200,000 to 3,000,000 words and span roughly a factor of 
fifty in computational requirements. The forecast models 


use two basic data flow structures and are grouped as 
follows: the five and ten layer, coarse horizontal grid 

models, PECMCV and PECHFV, and the five and ten layer fine 
horizontal grid models, PEFHCV and PEFHFV. 

During the eoxirse of developing these forecast models, 
three different computer systems were utilized: a CDC 6500, 

a CDC CYBER 175 and a CDC CYBER 76, all using CDC software. 
Thus, an attempt at portability was required and largely 
realized, although the forecast models are in no sense 
totally portable between all computer systems with FOR- 
TRAN IV as their high level language. 

It is not the intent of this section to give complete 
detail down to the definition of the last local variable in 
each subprogram, but to describe the data structure of each 
program and to describe the code structure v;ith sufficient 
detail to be easily understandable without overwhelming rhe 
reader with trivia. Also, since there is a large amount of 
duplication of code between the four forecast models, dif- 
ferentiation between programs in the explanations which 
follow is only made when necessary. 

Section A explains the data structures of the four 
forecast models. Section B describes the program structures 
of the models, giving block diagrams, defining the various 
common blocks and giving a functional description of each 
of the subprograms. 
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A. Data Structure 

One of the most critical factors involved in attempting 
to design an efficient computer program is an efficient 
data storage scheme, since a primitive equation forecast 
model (PET4) requires the handling of large amounts of data. 
Section 1 gives a general picture of the data storage 
philosophy used in the models. Section 2 then discusses 
the details of the two versions of data storage structure 
used in the models giving diagrams, tables and illustrative 
examples of the working of this method of data storage. 

1, Data Storage - Philosophy 

Most PE models numerically approximate the solution 
of a complicated set of partial differential equations 
the. primitive equations — through the use of finite differ- 
ences. While these PEM are three-dimensional in space, the 
treatment of the third space dimension, the vertical dimen- 
sion, is quite different from the two horizontal dimensions. 
The use of the hydrostatic assumption to render the numerical 
solution feasible and vertical column computations in the 
heating packages immediately offer themselves as obvious 
examples of the difference between the horizontal and the 
vertical dimensions. 

If one assumes some sort of grid system, for 
example, either a polar stereographic projection or a 
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spherical polar grid, one then has an array of data, say MxN 
in the horizontal with K vertical levels. Since there v;ill 
be at least five variables at every point, the total number 
of pieces of information can be quite large, on the order 
of millions of words. Except for the very largest computer 
memories, this vast amount or data, together with the nec- 
essary computer code, cannot be conta ’ ned in memory. Thus, 
one is faced with the problem of de Ting some sort of 
efficient storage scheme for the data. 

Meteorologists are accust'iiied to looking at weather 
data on horizontal planes; for exfflnple, the geopotential 
heights at 500 mb. Very few outputs from a PIM are pre- 
sented as vertical slices. This, of course, is physically 
quite reasonable, since the scales of motion are quite 
different between the horizontal and the vertical, as are 
the grid distances in a PEM, Hov/ever, horizontal planes 
are not the best way to organize the data storage for com- 
putation in a PIM. 

Since it has been assumed that the numerical method 
used to approximate the solution of the primitive equations 
is finite differencing, it is easily seen that an entire 
horizontal plane is not necessarily needed in core simul- 
taneously. The second order finite difference approximations 
for horizontal derivatives used for the forecast models 
are given by 
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’f one imagines the j subscript being held fixed and i 
varying, then there are only three j rows involved in 
the computation along a single j row, the one being com- 
puted, the one below and the one above. In addition, since 
the heating package computations are all performed on ver- 
tical columns, it would seem reasonable to use a data 
storage scheme such as the one shown in Figure III-l. 


The scheme shown in Figure III-l has the computa- 
tional data broken down into vertical slices which are MxK. 
There are N of these slices. It has been seen that only 
three of these slices are needed at any one time to compute 
along a given j row. Thus, if one were to construct a file 
with N records on it, each containing one vertical slice, 
one could perform a forward integration of the primitive 
equations by sweeping this file in sequence. Of course, a 
new file should be constructed as this file is swept, in 
order to have data for the next step. The details of such 
a file scheme will be explained in the next section. 

A PEM consists of three basic sections, the ini- 
tialization section, the computational section which 


I 
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actually integrates the primitive equations, and an out- 
put section. By far the most computer time is spent in 
the computational section; thus, it seems reasonable to 
arrange the data storage for maximum convenience of this 
section of the code. This hypothetical PEM would take, 
as input, data on horizontal planes (pressure surfaces), 
since this is the way that analyses are performed. It 
would then interpolate this data to the sigma surfaces and 
finally sort out this data into vertical slices for use by 
the computational section. The computational section would 
then integrate the data forward in time to obtain a fore- 
cast and pass it to the output section. The output section 
will then sort out the data into horizontal planes and 
interpolate, smooth, etc, , for viewing by the users of the 
forecast. 

An obvious question that arises is, does this 
method of data storage result in an efficient PEM from the 
standpoint of central memory and computational speed? The 
answer is that it does and that this method of data storage 
is in use at numerous meteorological centers. As will be 
explained in the next section, rotating storage can be used 
quite effectively and more efficient utilization of central 
memory can be obtained when using the vertical slice storage 
technique. 
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2. Data Storage - Realization 



f 

t; 

i 

I 

I 


The concept of a data file consisting of N records , 
each representing a vertical slice of the three-dimensional 
crid system, was introduced in the previous section. If 
the grid is MxNxK, then there will be N records on the file. 
Each record will hold data for an MxK vertical slice. It 
was also pointed out that, in order to perform computations 
on the j row, only the j-1 and j+1 rows are needed in addi- 
tion to the j row. In this section, the file structure, 
memory storage requirements and data flow will be examined 
in detail. 

For the sake of clarity, let us call the existing 
data file lOLD and the data file created by one sweep of 
the data base INEW. If it is assumed that centered time 
differencing is used; i.e. , 

Ifl Ht [III-33 

and that the numerical technique of pressure gradient force 
averaging is used, then file lOLD will contain information 
at the n-1 and n time levels. After the first sweep of lOLD, 
INEW will contain partially updated information at the n and 
n+1 time levels (temperatures, moistures and geopotentials) 
along with the non-updated information. At the end of the 
sweep, the file INEW becomes the file lOLD and a second 
sweep of the data is made to integrate the momentum equations 


i 
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thereby completing the time step. The file INEW now con- 
tains all updated data at time levels n and n+1. This 
process is repeated until the desired forecast time is 
reached . 

In order to achieve more efficient use of central 
memory and to reduce the size of the data files, variables 
with time levels are ”time packed" on a 2:1 basis. That 
is, the n time level occupies the upper thirty bits of the 
word and the n-1 time level occupies the lower thirty bits 
of the word. Certain other infrequently referenced vari- 
ables are also packed on a 2:1 basis. A set of constructs 
which take the form of FORTRAN statement functions have 
been designed to be able to manipulate these packed data 
words. These statement functions are listed below with a 
brief explanation of their function: 

retrieving the n-1 time level 
or lower value - 

SNM(V) = SHIFT (V, 30) [III. 4] 

pack two values into a word - 

PACK(VNP,VN) = (VNP.AND.77. . .70. . .OB) .OR. 

(SHIFT (VN, -30) .AND.O. . .07. . .7B) [III. 5] 
store a value into the top half 
of a packed word - 

STOR(VNP,VN) = (VNP.AND.77. . .70. . .OB) .OR 

(VN.AND.O. . .07. . .7B) [III. 6] 
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Since this "time packing" method relies on the CDC 
FTN compiler intrinsic function SHIFT, the forecast models 
are not portable to other compilers. However, this was thought 
to be worthwhile since it affords a considerable compression 
in data file size and has a minimal effect on the computation 
time (approximately a 3% increase) . 
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a. Data Storage - PECHCV and PECHFV 

The concept of an old data file being read and a 
new data file being created on each sweep of the data base 
was presented in the previous section. In addition, the 
idea of "time packing" and the necessary methods of acces- 
sing the packed data were also presented. The next item 
to be considered is an efficient means of flowing the data 
to and from central memory in order to implement the com- 
putation . 

A diagram of the necessary central memory storage 
buffers for this data storage scheme to work is given in 
Figure III-2. Two buffer areas are required, AO for the 
three working rows and the row being input, and AN for the 
row being created and the row being output. It has been 
assumed that the computer hardware will allow simultaneous 
I/O and computation-buffered computation? thus, the row JIN 
in array AO and the row JM in array AN. Further, since it 
is time consuming to move data about in core, the buffers 
AO and AN are circular buffers, with the pointers redefined 
as the computation proceeds. 

One more construct is needed before a given piece 
of data can be accessed in the arrays AO and AN. This is 
a method of pointing to the particular variable that is 
required on the particular row that is desired. For com- 
putational convenience/ two subscript functions have been 
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FIGURE lll-2t CENTRAL MEMORY STOOGE 

BUFFERS REQUIRED FOR ROW 
STORAGE SCHEME USED IN 
MODELS PECHGV AND PECHFV 
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defined, one for data that exists on only one vertical 
level, and one for data that exists on numerous vertical, 
or K, levels. These subscript functions use the row 
pointers, pointers to the variables in a row which are 
known as offsets, and the vertical level, k, for the second 
type of subscript function. 

The two subscript functions, which take the form 

of FORTRAN statement functions, are defined below. 

ISUBS (J/IOFF) = J*m7AR+IQFF [III. 7] 

for data on only one vertical level and 

ISUBK(J,IOFF,K) = J*NVAR+IOFF+ (K-1) *IMAX [III. 8] 

for data on many vertical levels. In these definitions, 

J = Pointer to the row being referenced in the 
buffer array 

lOFF = Offset of the particular variable desired 
K = Vertical level desired 
NVAR = Total number of words per row record 
IMAX = Total length of a row (range on 
It is assumed that the pointers to the AO array ( JMN , JN, 

JPN and JIN) range between 0 and 3 and that the pointers 
to the AN array {JMNP and JNP) range between 0 and 1. These 
subscript functions then point to the 1=0 element for 
ISUBS and the 1=0 element on level K for ISUBK . 

A simple example of the use of the subscript func- 
tions is called for at this point. Assume that it is desired 
to compute on line J, level K; then the FORTRAN code for 
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doing this would be as follows: 

DO 10 K=1,KMAX 
IU=ISUBK(JN,UN,K) 

DO 10 I=2,IMAXM 

DO (I , K)= (AO (IU+I+1) -AO (IU+I-1) ) /TDX 
10 CONTINUE 

The raechanics for accessing data in the storage 
buffers AO and AN have now been constructed. The data access 
mechanism together with the buffer structure shown in Figure 
III-2 allows a buffered computation to be performed; that 
is, simultaneous computations and I/O to the files lOLD 
and INEW. The next, and final, step is to put these ideas 
together to generate the data flow for one complete time 
step. The data flow shown assiimes that free slip, insulated 
wall boundary conditions are used. Figure III-3 shows the 
data flow for one sweep of the data file for the forecast 
models PECHCV and PECHFV. 

Finally, Table III-l defines the offsets which are 
used by the two forecast models in the pointer scheme to 
access data in the buffer arrays AO and AN. 


Defxne Poxnters JMN, JN, JP 
for AO and JMNP. JNP for AN 


Brxng xn First 3 Rows from lOLD 
to AO, Lines JMN, JN 


Compute Lower Boundary Condxtxons 


opy Row JMN xn AO to Row JMNP 
N 


Start Buffer of Row from lOLD to 
AO, Row JIN, Except on Last Time 
Through 


Start Buffer of Row JMNP in AN 
to INEW 


compute Left-Rxght Boundary 
Conditions 


Compute Upper Boundary Conditions 
If Last Time Throuah 


Copy Row JN in AO to Row JNP in 
AN, Update Below 


ompute Diagnostics for Row JN 


Compute on Row JN in AO Storing 
Updated Data in Row JNP in AN 


mm 


Check Buffer of Row JIN xn AO from j 
lOLD, Except on Le'“t Time Through | 


Check Buffer of Row JMNP in AN to 
INEW 


FIGURE III- 3: DATA FI 7. 7 FOR ONE SWEEP OP 

THF DATA FILE - MODELS PECHCV 
ANl' PECHFV 
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TABLE III-l; OFFSETS OF VARIABLES IN ARRAYS 
AO AND AN FOR PECHCV AND PECHFV 


PiSCRIPTION 

Initial terrain pressure, ttq, 
and " reduction to sea level 
(packed) 

Map factor, m 

Geopotential height, 

U wind, n and n-1 time levels 

V wind, n and n-1 time levels 

Temperature, n and n-1 
time levels 

Moisture, n and n-1 
time levels 

Terrain pressure, n+] 
and n time levels 

Terrain pressure, u, n 
and n-1 time levels 

Terrain geopotential 

Alt e do and land- sea- ice 
indicator (packed) 

Sea surface temperature 

Sea surface vapor pressure 

Tendency truncator 

Sin of latitude, sin<(j 

Precipitation 


NAME 

5 

LAYER 

10 

LAYER 


Level No. 

Level 

No , 

PTIR 

1 

0 

1 

0 

MF 

1 

63 

1 

63 

PHI 

5 

126 

10 

126 

DN 

5 

441 

10 

756 

VN 

5 

756 

10 

1386 

TN 

5 

1071 

10 

2016 

QN 

3 

1386 

6 

2646 

PTNP 

1 

1575 

1 

3024 

PTN 

1 

1638 

1 

3087 

G2 

1 

1701 

1 

3150 

ALSI 

1 

1764 

1 

3213 

SST 

1 

1827 

1 

3276 

QSS 

1 

1890 

1 

3339 

RC 

1 

1953 

1 

3402 

SL 

1 

2016 

1 

3465 

PRE 

1 

2079 

1 

3528 


TABLE III-l: OFFSETS OF VARI/^LES (Continued) 


DESCRIPTION 

NAME 

5 LAYER 
Level No, 

10 LAYER 
Level No. 

Heating Rates 

HT 

3 

2142 

3 

3591 

Vertical velocities, w 

VN 

4 

2331 

9 

3780 

X direction pressure 
gradient force, n and 
n-1 time levels 

PXN 

5 

2583 

10 

4347 

y direction pressure 
gradient force, n and 
n-1 time levels 

pyN 

5 

2898 

10 

4977 


NVAR 

— 

3213 


5607 


NOTE; IMAX =63 


b. Data Storage - PEFHCV and PEFHFV 

The magnitude of the data file that must be proc- 
essed for the fine horizontal resolution forecast models 
is eonsiderably larger than for the coarse horizontal res- 
olution models. This rather severe].y restricts the type of 
data flow scheme that can be constructed, since the problem 
becomes one of fitting the computation into the machine (a 
CYBER 76 was used as the target machine) . 

In order to fit the computation in the machine, a 
split row was devised. A split row contains all of the 
variables for a vertical slice, but they are grouped into 
two segments. The lower portion of the split row contains 
the variables that are required on the three rows J-1, J and 
J+1 during computation on row J. The upper portion of the 
split row contains variables that are only needed on row J . 
This allows smaller central memory buffers since there need 
be only one full row and two partial rows in central memory 
s imul t aneous ly . 

A diagram of the necessary central memory and ECS 
or LGM storage buffers foi this scheme to work is given in 
Figure III-4. Thert, are three rows in central memory, a 
partial J-1 row, AJM, a full J row, AJ, and a partial J+1 
row, AJP. IN ECS or ECM there is a buffer array, AIN, for 
three full rows. This buffer array is a circular buffer and 
has pointers: JLM, JL and JLP pointing to rows J-1, J and J+1, 
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AIN array 

Pointers JLM, JL and JLP 


FIGURE III-4; STORAGE BUFFERS REQUIRED 
FOR ROW STORAGE SCHEME 
IN MODELS PEFHCV AND PEFHFV 
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respectively. Due to the way that the CDC system software 
works / access to the files lOLD and INEW is made through 
the central memory array AJ with a block transfer either 
to or from ECS or LCM also being required. 

One more construct is needed before a given piece 
of data can be accessed in the arrays AJM, AJ and AJP. 

This is a method of pointing to the particular variable 
that is required. The offset of the variable can be used 
to point to single vertical level data since there are no 
row pointers to the central memory arrays. However, for 
computational convenience, a subscript function has been 
defined for data that exists on more than one vertical level, 
or K level. This subscript function, which takes the form 
of a FORTRAN statement function, is defined below; 

ISUBK(IOFF,K) = I0FF+(K-1) *IMAX [III. 91 

where 

lOFF = Offset of the particular variable desired 

K = Vertical level desired 

IMAX = Total length of a row (range on I) . 

For single level data the offset of the variable points to 
the I = G element, and for multiple vertical level data the 
subscript function ISUBK points to the 1=0 element for 
level K. An example of the use of the subscript function 
was given in the previous section. 

The mechanics for accessing data in the central 
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memory arrays AJM# AJ and AJP have now been constructed. 

The data access mechanism together with the buffer structure 
shown in Figure III-4 allows the computation to sweep the 
data base on file lOLD while constructing an updated data 
base on file INEW. Figure III-5 combines these ideas and 
presents the data flow for one complete sweep of the data 
file in the forecast models PEFBGV and PEPHFV. The data 
flow shown assumes that free slip, insulated wall boundary 
conditions are used. 

Finally, Table III-2 defines the offsets which 
are used by the two forecast models to access data in the 
buffer arrays AJM, AJ and AJP. 
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Define pointers JLM, JL and JLP 
to array AIN in LCM 


Buffer in first 3 rows to AJ from 
lOLD, block transfer to rows JLM, 
JL and JLP of AIN in LCM 


Block transfer first 3 rows from 
AIN ia LCM to AJM, AJ and AJP 
in SCM 


Compute lower boundary 
conditions 


Block transfer AJM in SCM to row 
JLM of AIN in LCM 


Block transfer row JLM of AIN in 
LCM to AJ in SCM 


Buffer out first row from AJ to 
INEW 


Initialize line counter LINE=2 


Block transfer rows JLM, JL and 
JLP of AIN in LCM to AJM, AJ and 
AJP in SCM 


Compute E-W boundary conditions 


Compute upper boundary conditions 


Compute on row J in AJ 

































TABLE II I- 2: OFFSETS OP VARIABLES IN ARRAYS 

AJM, AJ AND AJP FOR PEPHCV AND PEFHPV 


DESCRIPTION 

N^E 

5 LAYER 
LeveJl No, 

10 

Level 

LAYER 

No« 

Initial terrain pressure « 
TTo/ and IT reduction to 
sea level (packed) 

PTIR 

1 

0 

1 

0 

Map factor# m 

MP 

1 

187 

1 

187 

Geopotential heights , $ 

PHI 

5 

374 

10 

374 

U wind# n and n-1 time levels 

UN 

5 

1309 

10 

2244 

V wind# n and n-1 time levels 

VN 

5 

2244 

10 

4114 

Temperature# n and n-1 
time levels 

TN 

5 

3179 

10 

5984 

Moisture# n and n-1 
time levels 

QN 

3 

4114 

6 

7854 

Terrain pressure# it, n+1 
and n time levels 

PTNP 

1 

4675 

1 

8976 

Terrain pressure# it# 
n and n-1 time levels 

PTN 

1 

4862 

1 

9163 

Terrain geopotential 

GZ 

1 

5049 

1 

9350 

Albedo and land- sea- ice 
indicator (packed) 

ALSl 

1 

5236 

1 

9537 

Sea sui^face temperature 
and sea surface vapor 
pressure (packed) 

SSTQSS 

1 

5423 

1 

9724 

Tendency truncator 

RC 

I 

5610 

1 

9911 

Sin. of latitude# sin<^ 

SL 

1 

5797 

1 

10 #09 8 

Precipitation 

PRE 

1 

5984 

1 

10,285 
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TABLE OFFSETS OF VARIABLES (Coatinued) 


DESCRIPTION 

NAME 

5 LAYER 
Level No. 

10 LAYER 
Layer No , 

Heating Rates 

RT 

3 

6171 

3 

10,472 

Vertical velocities, w 

W 

4 

6732 

9 

11,033 

X direction pressure 
gradient force, n and 
n-1 time levels 

PXN 

5 

7480 

10 

12,716 

Y direction pressure 
gradient force, n and 
n-1 time levels 

PYN 

5 

8415 

10 

14,586 


NVARS == 


5049 


9,350 


NVARL = 


9350 


16,456 


NOTE: IMAX = 187 





B. Code Structure 


This section describes the program structure of the 
forecast models PECHCV, PECHEV, PEFHCV and PEFHFV. It also 
defines the various files used by the programs, defines the 
common blocks that are used and describes each of the sub- 
programs that constitute the models. The models are all 
fairly similar and are quite similar by pairs; that is, the 
coarse horizontal grid five and ten level models and the 
fine horizontal grid five and ten level models can be grouped. 
Thus, an attempt has been made to describe all of the codes 
simultaneously. Of course, some differences are noted for 
certain routines and special purpose routines for a partic- 
ular group and are described separately. 

The codes are written in CDC FORTRAN, which is a dia- 
lect of FORTRAN IV, for the GDC FORTRAN Extended Compiler 
(0PT==2 level) . The two groups of programs have been designed 
to use rotating mass storage communicating with central 
memory through two I/O channels in the case of PECHCV and 
PECHFV, and rotating mass storage along with central memory 
and either ECS or LCM for the codes PEFHCV and PEFHFV. 


The forecast models are all fairly large codes. Since 
the various sections of the codes — initialization, inte- 
gration and output — are each fairly large, the codes have 
been overlaid to more efficiently use central memory. The 


sections that follow describe the main divisions of the 
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forecast models by overlays. Section 1 describes the 
control overlay and gives a block diagram of the overlay 
structure, defines the I/O file structure and defines the 
common blocks that are global to all overlays. 

Section 2 describes the initialization overlay defining 
the common blocks local to this overlay and describing each 
of the subprograms that are used by this overlay. Similarly, 
Sections 3 and 4 describe the integration and output over- 
lays, respectively. 

1. Main Overlay - Programs PEGHCV, PECHPV, 

PEFHGV and PEFHFV 

The main overlay of the forecast model programs is 
very simple, since it only performs the function of declar- 
ing three global common blocks, defining various control 
variables and calling the actual working overlays. The 
overlay structure of the programs is shown in Figure III-6. 
Only one level of overlay is used for programs PECMCV and 
PECHPV, while two levels of overlay are required for the 
initialization and output sections of programs PEFHGV and 
PEFHFV due to their greater complexity. 

The various files used by the programs arc defined in 
Table III-3. Tables III-4 and III-5 define the variables 
that appear in the global common blocks CKTRL, and INDEX, 
respectively. Common block BUFPOINT consists of the offsets 
for the variables and contains the integer variables listed 







TABLE 

OUTPUT 

TAPEl 

TAPE2 

TAPE 3 

TAPE 4 

TAPES 

TAPES 


II I- 3: file STRUCTUBE OF 

FORECAST MODELS 


Used to write various diagnostic 
messages 

Input data file for initialization 
overlay; generited by analysis 
programs 

Random access file used by 
initialization and output 
overlays 

Restart data tape; contains a 
header record and contents of 
current lOLD file 

Output data file; contains fore- 
cast fields to be selectively 
plotted by graphics program and for 
use by the analysis programs 

(lOLD, INEW) integration data 

(lOLD, INEW) integration data 


TABLE III-4: DEFINITION OF VARIABLES 

IN COMMON BLOCK /CNTRL/ 


lOLD 

INEVV 

ISW 

ID (2) 

LIMI 

LIM 

ITAU 

ITSPHR 

NN 

NNOUT 

LCO 

L12 

DAY 

HR 

DT 

TDT 

MATS 


File name of integration data 
(TAPES or TAPE9) 


File name of integration data 
(TAPES or TAPE9) 


=1, SSWl is "on", signifying normal 
start with data from TAPEl 
=2, SSWl is "off", signifying restart 
with data from TAPE3 

1st word, forecast identifier (date-time group) 
2nd word, forecast data field identifier 

Starting cycle 

Number of integration and output cycles 
(2 X number of days + 1) 


Forecast time indicator (hours) 

= 0 during initialization and first output cycle 
= 12*LC0 for normal integration cycle 

Number of integration time steps per hour 

Time step counter, reset after each output 
cycle, also NN=0 signifies a Matsuno time step 

Number of time steps between normal 
output cycles, usually equal to 
12 X ITSPHR 


Output cycle counter, updated by 1 each 
output cycle 

Indicator set at 12-hour output cycle, 
used to zero rainfall 

Number of days from January 1 
of the forecast 

Hour of the day, GMT 

Time step, At 

Twice the time step, 2 At 

Logical variable used during Matsuno time step 
= .FALSE, during first half 
= .TRUE, during second half 
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TABLE III-5: DEFINITION OF VARIABLES 

IN COMMON BLOCK /INDEX/ 

ICL Start of computation on data row (=2) 

ICU End of computation on data row 

(=*IMAX-1) 

IMAXM Length of data row -1 

IMAX Length of data row 

JMAXM Number of data rows -1 

JMAX Number of data rows 

KMAXQ Number of moisture levels 

KMAXM Number of vertical levels -1 

KMAX Number of vertical levels 

LINE Number of row being worked on 

NFLD Size of a data record on TAPE2 (63 x 63) 

NVAR* Size of a data record on TARE 8 and TARE9 

* The above table is correct for codes PECHCV 
and RECHFV. For codes PEFHCV and REFHFV the 
last line is replaced by 

NFLDEX Expanded data record (187 x 187) 

NVARS Size of the AJM and AJP buffers 

NVARL Size of the AJ buffer and a data record 
on TAPES or TAPE9 
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in Table III-l for PECMCV and PECHFV, and Table III-2 
for PEFHCV and PEFHFV. 

2. Initialization Overlay Structure 

The function of the initialization overlay is to 
take as input data the fields generated by the analysis 
programs and to generate an initial state for the inte- 
gration overlay which actually performs the forecast. That 
is, the interpolation and synthesis of data fields to the 
PEM sigma surfaces using data produced by the analysib 
programs as described in Section II-5 are performed by 
the initialization overlay. 

Figure III-7 shows the structure of the initiali- 
zation overlay, giving the names of all the subprograms 
which constitute this overlay. 

Computer core storage requirements are minimized for 
the 187 X 187 grid model by dividing the initialization 
overlay into a primary and three secondary overlays . 

Figure III- 8 shows the overlay structure and the sub- 
programs included in the overlay. With the exception of 
a simple executive, INITEX, and some interpolation rou- 
tines, these subprograms are identical in both function and 
name to those of the 63 x 63 grid model. The following 
sectiors describe these routines and point out differences 
required for the two overlay structures. 
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FIGURE III-8 


INITIALIZATION O^TERLAY STRUCTOiRE 
FOR PROGRAMS PEFHCV AND PEFHFV 













The Icibeled common block , INCON, is local to this 
overlay and its variables are described in Table III-6 . 

Table III-7 defines the data fields on TAPEl which 
are produced by the analysis programs. Table III-8 defines 
the fields on TAPE2 which are produced both by the initiali- 
zation process to define the initial state of the forecast 
and by the output overlay to present the results of the 
forecast. 

a. Program INITEX 

Program INITEX is the primary overlay for the 
programs PEFHCV and PEFHFV. It performs a simple executive 
fxmction by calling overlay RSTART to initialize from a 
previously written data tape; or by calling the overlays 
IRITPE and SORTIT to perform initialization using the anal- 
ysis data. 

b. Program INITPE 

INITPE exists as the main program in the initiali- 
zation overlay for programs PECHCV and PECHFV. It exists 
as a secondary overlay for programs PEFHCV and PEFHFV. The 
program extracts data generated by the analysis programs 
from the sequential file TAPEl, Using these data, it xnter- 
polates and synthesizes the data fields on the sigma sur- 
faces and writes these new fields to the random file TAPE2. 
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TABLE III-6: DEFINITION OF VARIABLES 

IN COMMON BLOCK /INCON/ 


PLVLU3) 

Pressure levels of input analysis data: 
1000, 950, 900, 850, 700, 500, 400, 300, 
250, 200, 150, 100, 50 

PLVLN(13) 

fn {PLVL(K)1, K = 1,KMAX 

SIG(5)or(10) 

Sigma levels 0.9, 0.7, 0.5, 0.3, 0.1 
for PEFHCV and .95, .85, .75, .65, .55, 
.45, .35, .25, .15, .05 for PEFHFV 

ZLN(5)or (10) 

An [SIG(K)], K = 1,KMAX 

SLN(5)or (10) 
6 

R*An(©.9) for K-1, R*An[SIG(K)/SIG(K-l) ] ; K=2,KMAX 

2 

Gravitational acceleration = 9.80616 m/sec 

GRDC 

Constant used in calculation of sin (latitude) = 
(31.205)^ 

CTOK 

Conversion constant relating degrees 
Celsius to degrees Kelvin = 273.16 

AE 

Constant used in calculation of saturation 
vapor pressure = 21.656 

BE 

Constant used in calculation of saturation 
vapor pressure = 5418 
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TABLE III-7: CONTENTS OP THE INPUT 

DATA PILE - TAPEl 
C?E?\TED BY THE ANALYSIS PROGRAMS 


Data Field 


Geopotential heights of the 1,000 mb surface 


Temperatures on the 1,000 mb surface 

II n 95Q II 


Land>sea-ice indicator 
Surface pressure, Pg 
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TABLE III-7; INPUT DATA PILE - TAPEl (Continued) 


Record 


Data Field 


29 

30 

31 

Terrain height 

Sea surface temperature 

Albedo 

32 

U 

wind 

on the 1,000 mb 

surface 

33 

V 

II 

1,000 

II 

34 

U 

tl 

950 

II 

3S 

V 

11 

950 

II 

36 

U 

II 

900 

11 

37 

V 

II 

900 

11 

38 

0 

II 

850 

II 

39 

V 

It 

850 

II 

40 

U 

11 

700 

tl 

41 

V 

11 

700 

II 

42 

U 

II 

500 

II 

43 

V 

tf 

500 

II 

44 

U 

11 

400 

11 

45 

V 

11 

400 

11 

46 

U 

11 

300 

11 

47 

V 

It 

300 

11 

48 

U 

tl 

250 

It 

49 

V 

II 

250 

II 

50 

U 

11 

200 

II 

51 

V 

II 

200 

M 

52 

U 

II 

150 

II 

53 

V 

tl 

150 

ft 

54 

U 

» 

ito 

II 

55 

V 

O 

100 

11 


NOTE: Geopotentlal heights are in meters; 

temperatures are in degrees centigrade; 
surface pressure is in mb; 
terrain height is in meters; 
winds are in meters per second* 

All fields are 63 x 63 arrays. 
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TABLE 


CONTENTS OF RANDOM ACCESS 
PILE « TAPE2 

10 SIGMA LE\^L VERSION OF 
THE FORECAST MODEL 


Data Field 


heights of 1000 mb surface 
950 

. 9050 : : “ ^ 

om ^ 

. 700 

'■.-v-^-sbo;.-. 

400 ■■ ' 

200 

150 

50 


of 1000 mb surface ( “ C) 
: 950 
900 
850 
700 
500 
400 
• 300 
250 
200 
150 

i6b 

50 • 


Analyzed svucface pressure,, Fs (mb). 


Terrain height, 


(meters) 


Sea surface . temperature . ; I 
Land-sea-ice indicator 
Sin of latitude, sine 
Map factPf , m 
Tendency truncatof = 

Albedo. 

Terrain geopotential , gZiji 
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TABLE III-8: RANDOM ACCESS FILE - TAPE2 

{Continued) 


Record Data Field 


41 

Temperature on a = 

0.95 surface 

("K) 

42 

11 

0.85 


43 

It 

0.75 


44 

IT 

0.65 


45 

II 

0.55 ” 


46 

11 

0.45 


47 

11 

0.35 


48 

II 

0.25 


49 

II 

0.15 " 


50 

IT 

0.05 


51 

Vapor pressures on 

a = 0.95 surface (mb) 

52 

II 

0.85 


53 

II 

0.75 


54 

II 

0.65 


55 

II 

0.55 


56 

II 

0.45 '• 


57 

Terrain pressure, tt 

(mb) 


58 

59 

60 

Pressure difference 
in (tt ) 

Forecast Ps (mb) 

, Ps-TT (mb) 


61 

Geopotentials of e 

=* 0.95 surface 

(g« meters 

62 

11 

0.85 


63 

It 

0.75 


64 

II 

0.65 '• 


55 

II 

0.55 


66 

IT 

0.45 


67 

It 

0.35 


68 

II 

0.25 '• 


69 

11 

0.15 ” 


70 

II 

0.05 ” 


71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

U winds on 90''' mb surface 
" 700 

•• 500 

H 300 

" 250 

'• 150 

» 100 

(m/sec) 
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TABLE II I- 8: 


RANDOM ACCESS PILE 
(Continued) 


TAPE 2 


Data Field 

V winds OB 900 mb surface 
•' 700 

" 500 

" 300 

" 250 

'• 150 

•' 100 


(m/sec) 


U winds on ex 


(m/sec) 


V winds on cr = 


(m/sec) 


Forecast geopoten 


0.95 surface 
0.85 

0.75 " 

0.65 

0.55 

0.45 

0.35 

0.25 

0.15 

0.05 

0.95 surface 

0.85 

0.75 

0.65 

0.55 

0.45 

0.35 

C.25 

0.15 

0.05 


tial height of 1,000 mb surface (meters) 
" 950 

" 900 " 

" 850 " 

" 700 

■> 500 " 

" 400 " 

" 300 " 

.1 250 

" 200 ” 

” 150 •' 

•• 100 " 


III-42 
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Record 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 


TABLE III-8: RANDOM ACCESS FILE - TAPE2 

(Continued) 


Data Field 


Forecast temperatures on 1,000 mb surface (®C) 


II 


950 

tl 


If 


900 

»' 


II 


850 

II 


I! 


700 

11 


IT 


500 



11 


400 

II 


If 


300 

M 


It 


250 

11 


It 


200 

11 


II 


150 

II 


It 


100 

II 


Forecast 

U winds on 1,000 

mb surface 

(ro/sec) 


11 

950 

11 



tl 

900 

If 



II 

850 

II 



If 

700 

IP 



II 

500 

II 


: 

II 

TT 

400 

n n 

It 

11 


i 

11 

:>U U 

250 

It 


5 

II 

200 

IT 


j 

] 

II 

150 

M 



II 

100 

II 



Forecast 

V winds on 1,000 

rob surface 

(ro/sec) 

t 

II 

950 

II 


•i 

i 

II 

900 

It 



If 

850 

M 


1 

II 

700 

Tl 


II 

500 

II 


! 

j 

It 

400 

II 


11 

300 

II 


3 

II 

250 

II 


] 

II 

200 

tl 


i 

i 

II 

150 

II 


IT 

100 

It 


\ 

1 

1 

s 

1 

Sea surface saturation vapor pressure 

(mb) 

Forecast precipitation 

(cm) 


.1 

Vapor pressure on 1,000 mb surface 

(mb) 

1 

It 

900 

II 


tl 

700 

11 


3 

11 

500 

11 


1 

It 

300 

II 


\ 
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c. Subroutine SETPHIS 

Subroutine SETPHIS computes the initial geopoten- 
tials of the sigma surfaces. Input to the routine consists 
of the temperatures on the sigma surfaces and the moistures 
on the lower sigma surfaces. The geopotentials are computed 
using the form of the hydrostatic equation given by Equa- 
tions [11.107] and [11.108] given in Section II-D-1. 

Finally, the computed geopotentials are stored on the 
random access file TAPE2 for use in the forecast. 

d. Subroutine MOIST 

Subroutine MOIST computes the relative humidity 
on a pressure surface using the procedure based on the rel- 
ative vorticity as described in Section II-E-6. Input to 
the routine consists of the geopotential heights of the 
pressure surface. The routine then computes the relative 
vorticity on the pressure surface and, using the empirical 
relation given by Equation [11.133], computes the relative 
humidity on the pressure surface. 

Program INITPE calls this subroutine four times 
to generate moistures on the 1,000, 900, 700, 500 and 
300 mb surfaces. Data on these pressure surfaces is then 
used to obtain initial moisture data on the lower sigma 
surfaces through logarithmic pressure interpolation. 

These data fields are stored on the random access file 
TAPE2 for use in the forecast. 
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e. Subroutine TREAD 

Subroutine TREAD acts as a generalized input data 
routine to obtain analysis data from TAPEl. A storage array 
address and the data field ID are input to TREAD from the 
cr.lling routine. TREAD then searches a table (array ITP) 
to find the location of the desired data field and reads 
the data field from TAPEl into the storage array. Various 
data field orderings on TAPE! may be accommodated simply 
by changing the information stored in array ITP. 

f. Subroutine DATAIO 

Subroutine DATAIO is a generalized routine for 
writing and reading data fields to and from the random 
access file TAPEr:,, Two entry points exist for this routine: 
DATAOOT, which uses the system routine WRITHS for writing to 
the file, and DATAIN, which uses the system routine README 
for reading from the file. 

g. Subroutine RSTART 

In PECHCV and PECHEV, RSTART is a subroutine called 
by INITPE. In PEFHCV and PEFBFV, RSTART is a secondary 
overlay called by INITEX. Subroutine RSTART reads the re- 
start data tape, TAPE3, when a forecast run is being re- 
started and initializes the forecast model in order to 
continue the forecast. Restarts are controlled by sense 
switch 1. If SSWl is ”0N" , then a normal run is being made; 
if SSWl is "OFF", then a restart run is being made. 
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A restart data tape is written during each output 
cycle so, if the computer goes down or if for some reason 
a forecast is to be continued at a later time, the forecast 
can be resumed from the point at which the last restart 
tape was written. Table III-9 gives the format of the 
restart data tape TAPES . 

Thus subroutine RSTART reads in the header record 
from TAPES and initializes the required control variables. 
Next the routine reads in the forecast data records and 
recreates file lOLD. The forecast run will now proceed from 
this point in the forecast. 

h. Subroutine SORTIT 

In PECHCV and PECBFV, SORTIT is a subroutine 
called by INITPE. In PEFHCV and PEFHFV, SORTIT is a sec- 
ondary overlay called by INITEX. Since the output of the 
analysis programs and t’ - initialization routines are data 
on horizontal planes and the integration overlay performs 
computations on data arranged in vertical slices, an inter- 
face routine to properly structure the data is required. 
Subroutine SORTIT takes the initial data in horizontal pleine 
form which is stored on the random access file TAPE2 and 
constructs the forecast data file lOLD in vertical slice 
form for use by the integration overlay. 
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TABLE I I I- 9; RESTART DATA TAPE 
(TAPES) FORMAT 

1) Header record (tv;enty words) - 

contains information to identify the forecast 


1 

ID(IO) 

2 

ITAD 

3 

lODT 

4 

LCO 

5 

L12 

6-10 

Unused 

11 

DAY 

12 

HR 

13-20 

Unused 


2) Forecast data records - 

these contain the contents of file lOLD at ITAD 
(JMAX records in vertical slice format) 
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Due to core storage limitations, more than one 
pass must be made through the file lOLD to completely ini- 
tialize the forecast. In PECHCV and PECHFV, each pass 
consists of reading a maximum of ten planes of data from 
TAPE2 and. then sweeping file lOLD from row one to row JMAX, 
writing out the initial data. In PEFHCV and PEFHFV, each 
pass consists of reading a maximum of twenty planes of data 
from TAPE2 to central memory and then block transferring 
the data to LCM temporarily. Next, the file TOLD is swept 
so that on each pass the data fields are block transferred 
back to central memory, expanded from a 63 x 63 grid to a 
187 X 187 grid and written to the file lOLD a row at a time. 
The expansion is accomplished using a bi-cubic spline 
interpolation . 


i. Interpolation Routines 

The interpolation routines are called by gOR'riT 
in programs PEFHCV and PEFHFV. The routines perform a 
bi-cubic spline interpolation to expand the data fields 
from a 63 X 63 to a 187 x 187 grid. The names of these 
routines are as follows; 


IBCIEU 

ICSEVU 

ICSICO 

UERTST 


These routines have been extracted from the International 


Mathematical and Statistical Library (IMSL) . 
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3. Integration Overlay - Program INTGPE 

The function of the integration overlay is to numer- 
ically integrate the primitive equations; that is, to 
step the data fields forward in time to produce a forecast. 
Input to this overlay consists initially of the data field 
specification produced by the initialization overlay and, 
subsequently, the current state of the forecast contained 
on file lOLD. After initialization, a forecast consists of 
cycling through the integration and output overlays until 
the desired forecast time has been reached. 

Figure III- 9 shows the structure of the integration 
overlay, giving the names of all the subprograms which 
constitute this overlay. The function of each of these 
routines is described in the sections that follow. 

Tables ill- 10 and III-ll define the variables that 
appear in the common blocks that are local to the inte- 
gration overlay for programs PECHCV - PECHFV and PEFHCV 
- PEFHFV, respectively, TABLE iII-12 defines the variables 
in common block /CON/ which is the same for all four fore- 
cast model programs. 
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TABLE III- lOj DEFINITION OF VARIABLES IN 
COMMON BLOCKS LOCAL TO THE 
PROGI'AMS PECHCV AND PECHFV 


Common 


Common 


Common 


Block /BUFF/ 


AO(NVAR,4) 


Circular input buffer for rows 
JMN, JN, JPN and JIN 


AN(NVAR,2) 


Circular output buffer for rows 
JMNP and JNP 


Block /TEND/ 

OT(63,KMAX) 


VT(63,KMAX) 


TT(63,KMAX) 


, tendencies in tru for all 
levels on row JN 


3ttv 

9t 


tendencies in 


levels on row JN 


TTv for all 


f tendencies in 
levels on row JN 


1TT for all 


QT(63,KliAXQ) 


PIT (6 3) 


, tendencies in irq for all 
levels on row JN 

3 tt 

, tendency in terrain pressure, , 
on row JN 


Block /LINEAR/ 

DPIDX{63,KMAX) Kurihara modification pressure 

gradient terms in x direction for all 
levels on row JN 

DPIDY (63,KMAX) Kurihara modification pressure 

gradient terms in y direction for all 
levels on row JN 

W(63,KMAXM) Vertical velocities for all levels 
of row JN 


Common 


Note 


TABLE III- 10: 


Block /LOCAL/ 
UPM(63) 


VPM(63, 3) 


Definition of Variables (Continued) 


— at a particular level on row JN 
m 


— at a particular leve"' on rows 
JMN, JN and JPN 


PILN(63,3) ln(?r) on rows JMN, JN and JPN 

PIKAP(6 3) on row JN 

TCOR(63 , 3 ,KMAX) temperatures corrected for vapor 

pressure at all levels on rows JMN, 
JN and JPN 


AB(63,KMAX) 

QC(63,KMAXQ) 

QCS(63,KMAXQ) 

TC(63,KMAX) 


NOTE : TCOP. is equivalenced to the 

following variables : 

(6a + 6g) at all levels on row JN 

vapor pressure, q, at all levels on 
row JN 

saturation vapor pressure, q^, at all 
levels on row JN 

temperature at all levels on row JN 


PC(63,KMAXQ) 

RF(63) 


pressure at moisture levels on 
row JN 

precipitation on row JN 


Symbolic dimensions are given where appropriate since 
the size of some arrays varies with the forecast model. 
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TABLE III- 11; DEFINITION OF VARIABLES IN 
COMMON BLOCKS LOCAL TO THE 
PROGRAMS PEFMCV AND PEFHFV 


Coirnnon Block /BUFFL/ 
AIN(NVARL,3) 


Circular input buffer for rows JM, 

J and JP; resides in either ECS or LCM 


Common Block /BUFFS/ 
AJM(NVARS) 


AJ (NVARL) 
AJP(NVARS) 


Common Block /LOCAL/ 
PIKAP(187) 
W(187,KMAXM) 


Buffer for partial row JM in central 
memory 

Buffer for full row J in central memory 

Buffer for partial row JP in central 
memory 


TT on row J 

Vertical velocities for all levels on 
row J 


DUM(187,KMAX+4) 


a) PILN(187,3) 

b) UPM{187) 


Dummy array 

NOTE: DUM is equivalenced to the 

following variables 

iln(ir) on rows JM, J and JP 

on row J at a particular level 


VPM(187,3) 


AB(187,KMAX) 
c) PIF(187,3) 
PEPS (187) 


* — on rows JM, J and JP at a particular 
m 

level 

2 

4cm (6a-t<S8) at all levels on row J 
2Ax 

TT reduction factor to sea level 

TT change due to nonlinear smoothing 
operator 
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TABLE III-ll; Definition of Variables (Continued) 


Common Block /LTEND/ 

TENDP (187, 2*KMAX) Dummy tendency array 


DPIDX(187,KMAX) 


Pressure gradient force terras in x 
direction for all levels on row J 


DPIDY (187,KMAX) 


Pressure gradient force terms in y 
direction for all levels on row J 


TCOR(187, 3,KMAX) 


Temperatures corrected for vapor 
pressure at all levels on rov?s JM, 
J and JP 


a) PIT (187) 


TT(187,KMAX) 


NOTE : /LTEND/ resides in central memory 

on the Cyber 175 and in LCM on the 
Cyber 76 

NOTE ; TENDP is equivalenced to the 
following variables: 

3tt 

, terrain pressure tendency on rov; J 


9ttT 
9t ' 
row J 


tendencies in ttT for all levels on 


QT (187,KMAXQ) 


9t ' 

on row J 


, tendencies in irq for moisture levels 


b) DT(187,KMAX) 


9t ' 
row J 


, tendencies in ttu for all levels on 


VT (187,KMAX) 


9t ' 
row J 


, tendencies in ttv for all levels on 


NOTE i TCOR is equivalenced to the 
following variables: 


PC (187,KMAXQ) 


Pressure at moisture levels on row J 
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TABLE III-ll: Definition of Variables (Continued) 


QC(187,KMAXQ) 


QCS(187,KMAXQ) 


TC(187,KMAX) 


Vapor pressure at moisture levels on 
row J 

Saturation vapor pressure at moisture 
levels on row J 

Temperature at all levels on row J 


NOTE: Symbolic dimensions are given where appropriate since 

the size of the array varies with the forecast model. 
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TABLE I II- 12: 
D 

DSQ 

TD 

OFD 

DSIG 

TDSIG 

SIG(KMAX) 

CT (KMAXM) 
SLN (KMAX) 

S9* 

S7* 

S5* 

S3* 

SI* 

AE 

BE 

G 


DEFINITION OF VARIABLES IN COMMON BLOCK /CON/ 


Ax = Ay = 381000 for coarse mesh models and 
127000 meters for fine mesh models 

Ax^ = Ay^ 

2Ax = 2Ay 
1/4 Ax 

Ao = 0.2 for five layer models and G.l for 
ten layer models 

2A(j 

Si^a levels 0.9, 0.7, 0.5, 0.3, 0.1 for five 
layer models, and 0.95, 0.85, 0.75, 0.65, 0.5 
0.45, 0.35, 0.25, 0.15, 0.05 for ten layer 
models 

In t SIG{K)/SIG(K-1)3 ; K = 2, KMAX 

R-ln(0.9) or In (0.95) for K=l, R-CT{K) for 
K=2 , KMAX 

In (0.9) 

In (0.7) 

In (0.5) 

In (0.3) 

In (0.1) 

Constant used in calculation of saturation 
vapor pressure = 21.656 

Constant used in calculation of saturation 
vapor pressure = 5418.0 

2 

Gravitational acceleration = 9.80616 m/sec 
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TABLE III-12; Definition of Variables (Continued) 


R 

XKAP 


EC 


CD 

DK 


CK(KMAX) 


Gas constant = 287.04 
K = ^/Cp = 0.2858 

-4 

Coriolis parameter = 29. = 1.4584 2x10 
Drag coefficient = 0.0015 

6 2 2 

Diffusion coefficient = 1x10 m /sec for 

coarse mesh models and Ixio"^ for fine mesh 
models ^ ■■ 

[1000/SIG(K)] **XKAR,K=1,KMAX 


* These five constants are appropriate to the five layer 
models. They are replaced with S95=ln (0 . 9 5} , etc. 
for the ten layer models. 

_ ^ ; 

Note; Symbolic dimensions are given where appropriate 
since the size of the arrays varies with the forecast models. 


j 


r- 

s 

!■ 
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a. Subroutine DEFINE 

Subroutine DEFINE has the simple function of 
defining various constants that are used in the integration 
overlay and zeroing the buffer arrays of this overlay. The 
subroutine is called once per integration-output cycle upon 
entering the integration overlay. 

b. Subroutine LSTEPP 

Subroutine LSTEPP performs the first half of a 


pressure-gradient force averaged time step integrating 


forward one time step the pressure tendency, thermodynamic 
energy and moisture equations . The routine controls the 
data flow for one sweep of the data base and calls the 
routines necessary to compute diagnostics, boundary con- 
ditions, forcing functions and diabetic effects. Figure 
III- 9 shows the controlling function of subroutine LSTEPP. 

Figure III- 3 , which appears in Section III- 2-a 
shows the data flow for one sweep of the data base. This 
figure concentrates on the data flow aspect of the routine 
and the computational portion is denoted by the box labeled 
"computation on row J". Figure III- 10 is an expansion of 
this box for LSTEPP and shows the step-by-step computations 
and the subroutines that are called. 
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c. Subroutine PITEND 

Subroutine PITEND computes the Kurihara pressure 
gradient terms for use in the thermodynamic energy equation , 
the pressure tendencies and the vertical velocities. The 
terrain pressure tendency is computed using Equation [11.104] 
and the vertical velocity is computed using Equation [11,106]. 

d. Subroutine PISM 

Subroutine PISM performs the non-linear terrain 
pressure smoothing that is described in Section ll-D-2. The 
non-linear smoothing operator given by Equation [11.118] is 
applied along an entire row immediately after the terrain 
pressure has been advanced in time. 

e. Subroutine TQTEND 

Subroutine TQTEND computes the tendencies for the 
temperature and vapor pressure using Equations [11.105] and 
[11,109], respectively. 

The computation of the temperature tendency 
requires evaluation of the horizontal Arakawa advective 
operator, the vertical flux and the adiabatic heating/ 
cooling term. Computation of the vapor pressure tendency 
involves the horizontal Arakawa advective operator and the 
vertical flux. 


f. Subroutine LARCON 

Subroutine LARCON computes the amount of 
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precipitation resulting from the removal of excess moisture 
from super-saturated regions of the atmosphere and the 
temperature change due to the release of latent heat. The 
precipitation algorithm used in this computation is described 
in Section II-C-2. 

g . Subroutine HEATING 

Subroutine HEATING performs the diabatic heating 
computations which include shortwave and longwave radiation 
and sensible heating from the planetary boundary layer as 
described in Section II-C-1. Evaporation of moisture from 
the surface into the lowest sigma level is also computed. 

Since the radiation and boundary layer computations 
are complicated and time-consuming , they are only performed 
once every model hour (normally every four time steps , when 
a fifteen-minute time step is used) . Temperature tendencies 
are then computed and saved for use over the next hour of 
integration. 


h. Function ALON 

Function ALON computes the longitude of a Northern 
Hemispheric polar stereographic (I,J) grid point in radians 
ranging from 0 -2ir clockwise from the Greenwich Meridian. 

i. Subroutine CONVADJ 

Subroutine CONVADJ performs the moist and dry 
convective adjustment as described in Section II-C-3. The 


f 


moist convective adjustment consists of Arakawa’s param- 
eterized cumulus convection which may result in precipita- 
tion and adjustment of the temperature and moisture in the 
lower three sigma layers for the five layer models and in 
the lower six sigma layers for the ten layer models. The 
dry convective adjustment process consists of checking the 
model atmosphere for static stability and modifying the 
temperature structure when necessary to achieve stability. 

j . Subroutine SETi’HI 

Subroutine SETPHI computes the new geopotentials 
after all the prognostic variables have been updated using 
Equations tH*1071 and [11,108], 

k. Subroutine LSTEPV 

Subroutine LSTEPV performs the second half of a 
pressure-gradient force time averaged time step integrating the 
momentum equations forward one time step. The routine 
controls the data flow for one sweep of the data base and 
calls the routines necessary to compute diagnostics , boundary 
conditions and forcing functions. Figure III- 9 shows the 
controlling function of subroutine LSTEPV, 

Figure III- 3 , which appears in Section III- 2-a 
shows the data flow for one sweep of the data base. This 
figure concentrates on the data flow aspect of the routine 
and the computational portion is denoted by the box labeled 
"computation on row J". Figure III- 11 is an expansion of 
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this box for LSTEPV and shows the step-by-step computation 
and the subroutines that are called. 


1. Subroutine PGFTRM 


Subroutine PGFTRM computes the time averaged pressure 
gradient force terms as given by Equation [II. Ill] for use 
in the momentum equations. The Kurihara modification is 
used to evaluate the terms on local pressure surfaces rather 
than directly on the sigma surfaces. 

The routine also stores the vertical velocities 
into a working array for use in subsequent computations, 

m. Subroutine UVTEND 

Subroutine UVTEND computes the tendencies for the 
momentum equations using Equations [11.94] and [11.95]. 

The computation of the tendencies for the momentum 
equations involves the horizontal Arakawa advective operator, 
the vertical flux, the time averaged pressure-gradient force 
term, the Coriolis force and the diffusion term, 

n. Subroutine DRAG 

Subroutine DRAG modifies the u and v momentum ten- 
dencies in the lowest layer to account for frictional dissi- 
pation in the gross boundary layer. The friction terms are 
evaluated using Equations [11.98] and [11.99]. 
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o . Subroutine &NDRY 

Subroutine BNDRY computes the horizontal boundary 
conditions for the hemispheric polar stereographic grid. 

As explained in Section II-B, the horizontal boundaries con- 
sist of rigid impermeable vertical walls placed on the grid 
on the next to outermost row of grid points. The boundary 
conditions computed are given by Equation tll.9] . 

Three entry points exist for this subroutine: LBNDRY 

EWBNDRY and UBNDRY. LBNDRY computes values for the first 
row of the grid and is called once for each sweep of the 
data base. EVflBNDRY con.putes values for the end points of 
each row of the grid and is required for each row as the 
grid is swept. UBNDRY computes values for the last row of 
the grid and is called once for each sweep of the data base. 

p . Subroutine DIAGN 

Subroutine DIAGN computes and prints various layer 
mean diagnostic quantities in order to monitor the forecast. 
Table III-13 defines the diagnostic quantities computed by 
DIAGN. In addition to the mean values computed by the routine 
vertical column values are saved and printed for one arbi- 
trary point. 

There are three entry points in DIAGN: IDIAG, which 

initializes the diagnostic computations; DIAG, which computes 
diagnostics on a line-by-line basis; and DIAGP, which prints 
the layer mean diagnostics at the end of the time step. 
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TABLE III-13; DIAGNOSTIC QUANTITIES COMPUTED 
BY SUBROUTINE DIAGN 


i ^ 4 layer mean of the u wind 

A I J 


V, “ 2 V. , ,, layer mean of the v wind 


'k N i,j,k 

^ f J 


T, - i 2 T, 


i. - ^ ^ -1 V ' layer mean of the temperature 

^ # J 


^ ^ 'If' layer mean of the geopotential 

X I j 


~ h ^ * V' layer mean of the diabetic heating rate 

K 1^3, K 

X / J 

^ h ^ layer mean of the vertical velocity 


1/D 


M 2 vf mean of the vapor pressure 


1| J 


KE. = i IP... (u? . , + V? . .)/(2*R-T. . . ) , 

k N i ^ i»Dfk i/Utk i,T,k i,3/k 

X f J 


layer mean of kinetic energy 




layer mean of the square of the divergence 


VOUSQ^ - S j ^ '^i,j-i,k^^ 

X # J 


n - i 


layer mean of the square of the vorticity 
^ 2 ^^i j' hemispheric mean of the terrain pressure 

1 ' j ^ 


P = i £ 


^ 1 P. hemispheric mean of precipitation 

N X,D 


NOTE: N = (IC0-1CL+1) * (JMAX-2) 


ORIGlNAi; PAGE IS 
OF POOR QUAUry 
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q. Subroutine DIAGNL 

Subroutine DIAGNL computes and prints various 
zonal band diagnostic quantities in order to monitor the 
forecast. The zonal bands are ten degrees wide and range 
from the equator to the no:.th pole. Table III-14 defines 
the diagnostic quantities computed by DIAGNL. 

There are four entry points in DIAGNI<£ IDIAGL, 
which initializes the diagnostic computations ; DIAGLM, 
which computes all quantities except the eddy kinetic energy; 
DIAGLP, which computes the eddy kinetic energy; and DIAGPL, 
which prints the zonal band diagnostics at the end of the 
time step. 


r. Function POW 

Function POW computes the power function , to 

six-digit accuracy to achieve a savings in computer time 
and is used to compute tt^ and in the heating routine. 
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TABLE III-14; DIAGNOSTIC QUANTITIES COMPUTED BY 

SUBROUTINE DIAGNL 


If = ^ 2 ~ -i -i -i T' zonal mean of 

I<»K Nj^ l-»3 1»3»K J-fJ 

the tangential eoraponent of the wind 

V- „ = ^ 2 n. . vSinX, . + v. . ^^cosX. zonal mean of 

the normal component of the wind 

*®L,K ' RT 

It 1/3 

mean of the kinetic energy 


“I,K - It - Vk>' * 

It X tj 

zonal mean of the eddy kinetic energy 


i, j ,k 


)i 


^ zonal mean of the temperature 

Jj,K . . X,3,K 

It 1,3 

V? ^ L E |w] . . zonal mean of the vertical velocity 

L, K N . . . X, 3 

It 1,3 


DIVSQ. 


^ ^ t(u. 


L/K N- . “ . ^^“i+l,j,k - “i-l,j/k 

w D 

zonal mean of the square of the divergence 


^i/j+i/k ■ ''i/j-i/k^/^^''^ ' 


''°»sai.,K ■ r “''i+i,j,k - ''i-i.j.k " 

Ij 1 1 3 

zonal mean of the square of the vorticity 


NOTE: L denotes the zonal band and N_ denotes the number 

Lt 


of grid points lying in that band 


o> 
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4. Output Overlay - Program ODTPE 
The function of the output overlay is to take the 
results of the forecast contained on file lOLD and to 
produce meteorologically meaningful output; that is, to 
create surface pressure, precipitation and geopotential 
heights, temperatures and winds at standard levels for 
viewing by meteorologists. In addition, the output over- 
lay writes a restart data tape on each output cycle (nom- 
inally set at twelve-hour intervals). 

Figure III-12 shows the structure of the output over- 
lay, giving the names of all of the subprograms which 
constitute this overlay. Subroutine AVE is unigue to 
programs PEFHCV and PEFHFV; the remaining subroutines are 
functionally identical for all four of the prediction 
models. The function of each of these routines is described 
in the sections that follow. 

Table III-15 defines the fields which are produced 
by the output overlay and written to TAPE4. 

Table III-16 shows the common blocks local to this 


overlay 
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FIGURE 111-12: OUTPUT OVERLAY, PROGRAM OUTPE 



TABLE III-15: CONTENTS OP THE OUTPUT DATA FILE - 

TAPE4 PER OUTPUT CYCLE 


Record Data Field 

1 Forecast Surface Pressure, Ps (mb) 

2 Precipitation (cm) 

3 Geopotential Heights of 1,000 mb Surface (meters) 


4 

M 



II 

II 

950 

11 

It 


5 

II 



II 

II 

900 

II 

II 


6 

II 



II 

II 

850 

tl 

It 


7 

II 



II 

II 

700 

II 

II 


8 

II 



It 

It 

500 

It 

ir 


9 

II 



M 

IT 

400 

tr 

II 


10 

It 



tl 

II 

300 

tl 

II 


11 

IT 



M 

II 

250 

II 

11 


12 

II 



If 

11 

200 

11 

II 


13 

II 



11 

11 

150 

It 

11 


14 

11 



II 

II 

100 

II 

tl 


15 

Temperatures on 

the 

1,000 

mb 

Surface 


("C) 

16 

II 


II 

tl 

950 

II 

ir 


11 

17 

It 


IT 

Tl 

900 

11 

tl 


II 

T 8 

II 


11 

11 

850 

tl 

11 


It 

19 

11 


It 

II 

700 

II 

11 


11 

20 

II 


11 

If 

500 

11 

II 


tl 

21 

II 


II 

11 

400 

II 

11 


II 

22 

II 


11 

It 

300 

II 

II 


11 

23 

II 


IT 

It 

250 

11 

It 


II 

24 

tt 


tl 

II 

200 

II 

11 


II 

25 

II 


II 

tt 

150 

II 

II 


11 

26 

11 


It 

II 

100 

It 

II 


II 

27 

U winds 

on 

1000 

mb 

surface 

(m/sec) 



28 

V 

11 

1000 

Tl 

Tl 


II 



29 

U 

It 

950 

II 

It 


11 



30 

V 

II 

950 

II 

11 


It 
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TABLE III-15: 


OUTPUT DATA FILE 


TAPE4 (Continued) 


Record 

31 

U 

winds 

on 

900 

Data Field 
mb surface 

(m/sec) 

32 

V 

rr 

II 

900 

II 

II 

II 

33 

U 

n 

11 

850 

II 

II 

11 

34 

V 

M 

M 

850 

tl 

n 

It 

35 

U 

II 

II 

700 

II 

It 

II 

36 

V 

tl 

It 

700 

II 

II 

M 

37 

U 


It 

500 

M 

II 

tt 

38 

V 

n 

11 

500 

rt 

11 

It 

39 

U 

tl 

ft 

400 

tt 

11 

II 

40 

V 

tt 

II 

400 

II 

II 

II 

41 

u 

It 

11 

300 

II 

11 

11 

42 

V 

It 

II 

300 

tl 

11 

If 

43 

u 

11 

II 

250 

II 

It 

II 

44 

V 

M 

11 

250 

II 

II 

11 

45 

u 

M 

tl 

200 

11 

It 

It 

46 

V 

II 

II 

200 

M 

M 

II 

47 

u 

M 

II 

150 

II 

II 

II 

48 

V 

11 

11 

150 

11 

II 

ft 

49 

u 

II 

II 

100 

rt 

It 

II 

50 

V 

II 

M 

100 

11 

tt 

II 


NOTE: 


1 

1 

Each data field has an identification block of twenty | 

words appended to it, which uniquely identifies the I 

forecast base and forecast time plus parameters ] 

for plotting of the fields by the graphics package. | 
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TABLE III- 16: 


COMMON BLOCKS LOCAL TO 
THE OUTPUT OVERLAY 


Common blocks in programs PECUCV and PECHFV 


COMMON /WORK/FPl (63,63,10) 


COM^^ON/BUFFP/AO (3989) 


Used as temporary storage 
for the forecast fields on 
horizontal planes. 

Used as input/output buffer 
area for file lOLD. 


Common blocks in programs PEFHCV and PEFHEV 


COMMON/WORKX/DUMl (20) , 
FA(187,187) 


COMMON/BUFFP/AO (3969) 


Used as temporary storage for 
the forecast fields on hori- 
zontal planes and as an input/ 
output buffer area for file 
lOLD. 

Used as a working storage 
area . 
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a. Subroutine SOHTOOT 

Since the output of the forecast contained on 
file lOLD is data arranged in vertical slice format and 
the required output is data on horizontal planes, an inter 
face routine to properly structure the data is required. 
Subroutine SORTOUT takes forecast data in vertical slice 
format on sigma surfaces and constructs data on standard 
pressure surfaces and stores it on the random access file 
TAPE2 . 

In order to construct a data field for output 
purposes, the forecast values must be interpolated to 
standard pressure levels as the file lOLD is swept, rov? 
by row. This process creates horizontal plane data fields 
For graphical presentation, the surface pressure, geo- 
potential heights and temperatures are then smoothed to 
remove small-scale computational noise from the forecast 
values. The forecast fields are then written to the 
random access file TAPE2. Using subroutine PLTWRT which 
appends a twenty-word identification record to the fields, 
the forecast fields are written to the forecast file TAPE4 

b. Subroutine PLTWRT 

Subroutine PLTWRT appends a twenty-word identi- 
fication record to a data field and then writes the field 
to a specified output unit. The identification record 
contains forecast and field identifiers and contouring 
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information for the ^jraphical display package. 

e. Subromtine FILTR 

Subroutine FILTR is a non-linear smoothing routine 
used to remove computational noise from selected output 
forecast fields. Three types of filter can be chosen: a 

medium wave filter with a cutoff wave length of six grid 
lengths; a short wave filter v;ith a cutoff wave length of 
four grid lengths; and a long wave filter with a cutoff 
wave length of thirty grid lengths. 

d. Subroutine SMOOTH 

Subroutine SMOOTH is a diffusive type of data 
field smoother that only smooths points whose latitude 
is less than thirty degrees. These points are then smoothed 
according to the formula 

f! . = F. . + 0.1 (F) [III. 10] 

If”) 1 f j 

e. Subroutine PRT 

Subroutine PRT is a specialized printer contour 
routine which contours 63 x 63 data fields on four printer 
pages. It is used to output a selected few data fields 
to get a quick look at the quality of the forecast. 
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I 
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f . Subroutine RSTRTD 

Subroutine RSTRTD writes restart data to TAPE3 
during each output cycle for possible restart use in the 
event of computer malfunction or if it should bo desired 
to continue a forecast at a later date. The format of 
TAPE3 is given in Table III-9, which appears in Section 
III-B-2. 


g. Subroutine AVE 

Subroutine AVE appears in programs PEFHCV and 
PEFllFV. Its purpose is to transform forecast data from a 
187 X 187 grid to a 63 x 63 grid so that contours may be 
plotted using Subroutine PRT. Every third point along the 
boundaries of the 187 x 187 grid is moved to the boundaries 
of the 63 X 63 grid. The interior points are averaged via 
the following algorithm: 


1 ^ 
L 


F AVERAGE = , 4F. . + 2(F - . , . +F. +F. , . 

16 i 1,3 1+1,3 1»3+1 i“l»3 


^^i»j-l) ^ ^i+l,j-l "**^i+l,j + l ’’^^i-ljj+l 


+F 


i-1, j-1 


till. Ill 



i 

i 




■3 

i 
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